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'  This  is  the  final  report  in  a  series.  Overall  objectives  of  this  project 
are  concerned  with  developing  comprehensive  optical  and  thermal  signature  data 
bases,  the  development  and  evaluation  of  optical  and  thermal  canopy  radiation 
models,  and  the  interpretation  of  these  measurements.  Previous  technical  re¬ 
ports  in  this  series  have  described  optical  and  thermal  measurements  obtained 
over  a  coniferous  site  (Pinus  contorta)  in  Leadville,  Colorado.  This  earlier 
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'data  set  served  as  a  test  bed  for  the  development  and  initial  evaluation  of 
first,  individual  components,  that  is,  needle  and  leaves,  thermal  models,  and 
then  a  composite  canopy  terrain  model. 

The  objectives  of  the  work  reported  in  this  study  were  to  evaluate  the 
thermal  models  developed  under  a  wider  range  of  meteorological  conditions  and 
for  different  vegetation  types.  In  this  regard,  experiments  were  performed  on 
a  second  coniferous  site  (Pseudotsuga  manziesii)  near  Seattle,  Washington, -and 
a  deciduous  community  (oak-hickory)  at  the  Oak  Ridge  National  Laboratory, 
Tennessee.  As  part  of  the  evaluation  procedure  a  complete  sensitivity  analysi 
was  performed  for  the  model.  The  second  major  objective  of  the  study  reported 
here  was  a  restructuring  of  the  mathematical  model  which  enabled  a  factoring 
of  the  geometrical  characterization  of  the  canopy  in  terms  of  matrices  which 
can  be  convolved  with  the  energy  process  terms.  The  newly  structured  model 
more  easily  permits  the  precalculation  of  these  important  geometrical  char¬ 
acteristics  for  a  wide  variety  of  terrain  elements.  Finally,  two  parameter 
estimation  techniques  are  proposed  for  both  the  static,  steady-state,  thermal 
behavior  of  a  canopy  and  the  dynamic  or  time-dependent  implementation. 
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PREFACE 


The  research  described  in  this  report  was  conducted  by  personnel 
of  the  Department  of  Forest  and  Wood  Sciences,  College  of  Forestry  and 
Natural  Resources,  Colorado  State  University  (CSU)  from  1  October  1978  to 
1  February  1980  under  contract  No.  DACV!  39-77-C-0073  to  the  U.S.  Army 
Engineer  Waterways  Experiment  Station  (WES).  The  study  was  done  under 
Department  of  the  Army  Project  No.  4A762730AT42,  Task  A4,  Terrain/Oper¬ 
ations  Simulation,  Work  Unit  003,  Electromagnetic  Target  Surround  Char¬ 
acteristics  in  Natural  Terrains. 

Participating  project  personnel  concerned  with  the  tasks  described 
in  this  report  include  Dr.  James  A.  Smith,  Principal  Investigator;  Mr. 

K.  Jon  Ranson,  Research  Associate;  and  Mr.  Frank  Croft,  Graduate  Research 
Assistant.  In  addition,  very  significant  support  was  provided  by  Dr. 

Duong  Nguyen  of  the  Civil  Engineering  Department.  Dr.  Lee  Balick,  on 
assignment  at  the  WES  from  CSU,  was  responsible  for  the  technical  review 
of  the  report  and  numerous  suggestions  that  benefited  the  overall  quality 
of  this  report. 

Experimental  data  utilized  in  this  study  were  obtained  from  a  decid¬ 
uous  community  at  Oak  Ridge  National  Laboratory  in  conjunction  with  Dr. 

8.  Hutchison  of  the  Atmospheric  Turbulence  and  Diffusion  Laboratory  of 
the  National  Oceanic  and  Atmospheric  Administration.  Similarly,  measure¬ 
ments  were  obtained  over  a  Douglas-fir  community  in  cooperation  with  Dr. 

Leo  Fritschen  of  the  University  of  Washington.  Thermal  imagery  was  obtained 
by  the  Oregon  National  Guard  at  the  Washington  site. 

The  study  was  conducted  under  the  general  supervision  of  Dr.  John 
Harrison,  Chief  of  the  Environmental  Laboratory  (EL),  and  Mr.  Bob  Benn, 
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Chief  of  the  Environmental  Systems  Division,  EL.  Dr.  Lewis  E.  Link, 
Chief  of  the  Environmental  Constraints  Group,  EL,  was  Technical  Monitor 
for  the  study. 

Commanders  and  Directors  of  WES  during  the  conduct  of  this  study 
were  COL.  John  L.  Cannon,  CE,  and  COL.  Nelson  P.  Conover,  CE.  Technical 
Director  was  Mr.  Fred  R.  Brown. 


This  report  should  be  cited  as  follows: 

Smith,  J.A.,  Ranson,  K.J.,  Nguyen,  D.,  and  Link,  L.E.  1981. 
"Thermal  Vegetation  Canopy  Model  Studies,"  Technical  Report 
EL-81-6,  prepared  by  Colorado  State  University  in  collabora¬ 
tion  with  the  Environmental  Laboratory,  Waterways  Experiment 
Station,  for  the  U.S.  Army  Engineer  Waterways  Experiment  Sta¬ 
tion,  CE,  Vicksburg,  Miss. 
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THERMAL  VEGETATION  CANOPY  MODEL  STUDIES 


PART  I:  INTRODUCTION 

1.  This  technical  report  is  the  last  of  a  series  of  reports 
prepared  on  scene  radiation  dynamics.  Earlier  volumes  in  this  series 
have  described  the  development  of  models  for  optical  and  thermal  energy 
interactions  with  forest  and  grassland  vegetation  canopies.  Extensive 
field  measurement  efforts  done  in  cooperation  with  the  U.S.  Army  Engi¬ 
neer  Waterways  Experiment  Station  (WES)  have  also  been  separately  re¬ 
ported.  This  report  describes  further  efforts  in  thermal  model  devel¬ 
opment,  evaluation,  and  sensitivity  analysis.  Measurements  obtained 
over  a  Douglas-fir  (Pseudotsuga  menziesii)  experimental  site  near  Seattle, 
Washington,  and  an  oak-hickory,  deciduous  site  near  Oak  Ridge  National 
Laboratory,  Tennessee,  are  included.  At  both  sites  intensive  ground  instru¬ 
mentation  was  employed  as  well  as  thermal  overflights  provided  by  the  Oregon 
National  Guard  and  the  Georgia  National  Guard,  respectively.  In  addition, 
analyses  have  been  performed  with  data  from  Zweibriicken  Air  Force  Base  in 
the  Federal  Republic  of  Germany. 

2.  This  introduction  briefly  summarizes  the  following  topics  which 
are  explored  more  fully  in  the  body  of  the  report:  (a)  model  framework, 

(b)  sensitivity  analysis,  (c)  experimental  validation,  and  (d)  recom¬ 
mendations. 

Model  Framework 

3.  The  initial  thermal  canopy  model  utilized  in  this  study  is  de¬ 
scribed  in  the  report  by  Kimes,  Smith,  and  Ranson  (1979).  The  model  is  a 
plane-parallel  abstraction  of  a  vegetation  canopy  divided  into  three  hori¬ 
zontal  layers.  Furthermore,  steady-state  conditions  are  assumed. 
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4. 


An  energy-balance  formulation  of  the  model  may  be  given  in  vector 


/L 


form  by 


F(X,P,U)  =  0 


(1) 


where: 

F  =  (F1  F2  f3)  is  the  energy-balance  equation  for  layers  1,  2,  and  3, 
considering  the  following  energy  components:  longwave  transfers, 
shortwave  transfers,  sensible  heat,  and  evapotranspi ration 
X  =  (Xj  X2  X3)T  is  the  average  layer  temperature  vector  for  layers  1, 
2,  and  3 

P  =  (e.. ,  i=l,2,3  ,  i  =  l,2,3  eg  R1  A)  is  the  parameter  vector 

characterizing  the  canopy  layers 
e.j  ,cu  =  emissivity  and  absorptivity  of  the  vegetation  layer 
£g,c(g  =  emissivity  and  absorptivity  of  the  ground  layer 
R^  =  canopy  stomatal  resistance  to  water  vapor  diffusion 
S  =  longwave  flux  transfer  matrix  calculated  from  geometrical 
properties  of  the  canopy 

A  -  shortwave  flux  absorption  coefficient  vector 
U  =  (Tfl  Tg  WS  RH  SW)T  is  the  control  or  input  vector 
T  =  air  temperature 

a 

Tg  =  ground  temperature 
WS  =  wind  speed 
RH  =  relative  humidity 
SW  =  shortwave  flux 


5.  As  part  of  the  tasks  of  this  project,  £  was  rewritten  in  the 
following  form,  which  explicitly  factors  the  geometrical  properties  of  the 
canopy  from  the  remaining  energy  terms: 

F  =  *5  op  B(X)TS  +  B(X)  +  A  ■>  H(X)  +  LE(X)  (2) 


where: 

a  =  Stefan-Bol tzmann  constant 

£  =  vector  of  longwave  emission  terms 

H  =  vector  of  sensible  heat 

LE  =  vector  of  evapotranspiration  term 

The  significance  of  this  factorization  is  that  a  wide  variety  of  abstract 
or  canonical  canopies  may  be  characterized  by  precalculation  of  £  and  A 
matrices.  These  matrix  tables  may  then  be  convolved  with  the  appropriate 
meteorological  driving  variables  to  simulate  diurnal  behavior  for  a  wide 
spectrum  of  scenaries.  Five  standard  canopy  structures  of  three  different 
densities  are  given.  These  canopy  structure  combinations  represent  a  spec¬ 
trum  of  geometrical  structure-indexed  thermal  variations.  Other  combinations 
may  easily  be  calculated. 

6.  In  addition  a  view  factor  matrix  V£  is  precalculated  for  each 
canopy  characterization  which  is  used  to  calculate  thermal  exitance  W  as 
a  function  of  view  angle,  Q  . 

W(o)  =  VF(Layer,0)  BT  (3) 


where: 

W  =  the  predicted  canopy  exitance  at  view  angle,  0 


7.  Finally,  a  new  solution  of  the  energy-balance  equation  was 
formulated  utilizing  the  knowledge  of  the  £  function  which  permits  an 
explicit  evaluation  of  the  Jacobian. 

8.  Specifically,  a  modified  iterative  Newton-Raphson  technique  is 
employed  (Burden,  Faires,  and  Reynolds  1978). 

9.  Given  £  ,  U  for  a  given  time  period,  £(X,P,U)  becomes  a 
function  of  X  only.  Expanding  about  an  initial  guess,  ,  and  employing 
a  minimum  squared  error  criteria  yields 

■s* =  *  -  4  =  (j^rV  [-£(4)]  (4) 

where: 

J  =  the  Jacobian  evaluated  at  X  =  X  and  the  n+1  iteration  is 

—  -0 

given  by 


X 

-n+1 


X  +  <5X 

-n 


(5) 


Convergence  usually  occurs  within  a  few  iterations. 

10.  The  initial  guess  is  taken  to  be  air  temperature;  thus,  the 
solution  approach  may  be  interpreted  as  determining  the  modification  to 
the  air  temperature  profile  which  arises  when  a  canopy  is  inserted  into 
the  volume  space  under  consideration. 


Sensitivity  Analysis 

11.  A  sensitivity  analysis  was  performed  on  the  following  parameters 
and  input  variables: 

a.j  longwave  absorptivity  for  vegetation  layers  1,2,  and  3 

e .  longwave  emissivity  for  layers  1,2,  and  3 

£  ground  emissivity 
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R^  canopy  stomata!  resistance 

A.|  shortwave  absorption  in  vegetation  layers  1,  2,  and  3 
RH  relative  humidity 
Tg  ground  temperature 
WS  wind  speed 

T  air  temperature  above  the  canopy 

a 

Tac  a^r  temPerature  within  the  canopy 


Sensitivity  analysis  was  not  directly  performed  on  the  matrix  nor 
on  the  view  factor  matrix.  Rather,  the  above  analyses  were  repeated 
for  two  different  matrix  configurations.  One  corresponded  to  the 
Douglas-fir  canopy  and  the  second  to  an  oak-hickory  canopy. 

12.  Sensitivity  analysis  (Tomovic  1963)  involves  the  evaluation 
of  the  sensitivity  matrix: 


where : 


(6) 


=  layer  temperature  vector 
P  =  16-component  parameter/input  vector 


The  analysis  was  performed  in  each  case  for  x^,^  corresponding  to  a 
daytime  and  nighttime  representative  set  of  conditions. 

13.  The  first  order  perturbation  of  each  of  the  16  parameters  was 
evaluated  systematically,  solving  for  the  new  equilibrium  canopy  temperature 
profile  after  each  perturbation,  i.e., 


5X 


6P 


(7) 


t 


The  most  sensitive  parameter  of  the  model  was  found  to  be  the  air  temp¬ 
erature  within  the  canopy.  Next,  dependence  on  canopy  stomatal  resistance 
was  found  to  be  highly  nonlinear  for  the  low  values  of  .  The  depen¬ 
dence  of  canopy  temperature  on  most  other  parameters  was  found  to  be 
highly  linear. 


Experimental  Validation 

14.  Comparison  of  both  daytime  and  nighttime  measurements  for  the 
Doulgas-fir  and  oak-hickory  canopies  with  simulation  predictions  were 
carried  out.  For  both  of  the  canopies,  nighttime  simulations  deviated 
from  measured  values  by  2°C  or  less.  Daytime  simulations  underestimated 
measured  Douglas-fir  canopy  temperatures  by  a  maximum  of  2°C;  whereas, 
simulation  of  the  lower  canopy  for  oak-hickory  overestimated  temperatures 
by  a  maximum  of  4°C.  Deviation  patterns  could  be  explained  in  terms  of 
macroscopic  and  variable  environmental  conditions. 

Recommendations 

15.  Two  broad  categories  of  recommendations  are  made  in  the  enclosed 
report.  First,  several  suggestions  are  made  relative  to  improvements 
that  could  be  made  in  the  thermal  model  itself.  Secondly,  some  suggested 
approaches  for  estimating  required  parameters  in  the  model  from  observed 
data  are  given. 

16.  Sensitivity  analysis  has  indicated  the  importance  of  the  air 
temperature  within  the  canopy  as  an  input  to  the  model.  Further,  the 
validation  experiments  have  indicated  the  importance  of  utilizing  an 
appropriate  wind  speed  measurement.  Thus,  it  would  appear  to  be  appropriate 
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to  review  the  various  hypotheses  concerning  the  variation  in  air  temperature 
and  wind  speed  with  height.  The  model  is  easily  modified  to  include  a 
height  dependence  of  these  two  variables;  they  are  treated  as  constants 
simply  because  there  is  not  a  very  strong  rationale  for  choosing  among 
the  various  options.  In  a  similar  vein,  various  authors'  recommendations 
have  been  selected  for  analytic  representations  of  the  energy  budget  com¬ 
ponents.  It  may  be  useful  to  systematically  evaluate  several  alternative 
formulations.  Two  further  extensions  to  the  physics  of  the  model  would 
include  the  incorporation  of  a  ground  temperature  prediction  module  and  the 
expansion  of  the  steady-state  formulation  to  a  time-dependent  process,  that 
is,  allowing  for  heat  storage  within  the  canopy. 

17.  Finally,  further  analysis  of  the  structure  of  the  geometrical 
matrices,  that  is,  the  S  ,  A  ,  and  VF  matrices,  relative  to  the 
intrinsic  canopy  structure  variables  should  be  undertaken.  Specifically, 
the  possibility  of  further  factoring  these  matrices  in  terms  of  their 

leaf  area  index  dependence  and  their  dependence  upon  leaf  slope  distribution 
should  be  investigated.  It  may  be  possible  to  treat  the  density,  that  is 
the  leaf  area  index  dependence,  as  a  simple  scaling  influence  on  precal¬ 
culated  structural  forms.  If  an  analytic  decomposition  of  these  matrices 
in  terms  of  these  two  influences  is  not  possible,  numerical  approaches 
should  be  investigated.  A  faster,  more  tractable,  calculation  of  the 
shortwave  absorption  coefficient  should  be  given  high  priority. 

18.  Two  approaches  are  recommended  for  parameter  estimation  analysis. 
The  first  method  described  is  based  on  the  Kalman  filtering  techniques. 

The  linearization  of  the  model  in  terms  of  a  classic  state-space  framework 
is  outlined.  A  Kalman  filtering  approach  on  a  parameter  vector  or  an 
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augmented  state  vector  is  described  (Friedland  1972).  A  second  approach  to 
parameter  estimation  is  suggested,  which  is  based  on  the  use  of  sensitivity 
functions  (Durando  and  Leondes  1976).  This  approach  also  begins  with  a 
state-space  formulation  of  the  model  but  then  proceeds  to  use  the  sensitivity 
functions  to  calculate  an  unknown  parameter  vector  by  minimizing  the  square 
of  the  error  vector  between  predicted  and  measured  response. 

19.  The  appendixes  of  this  report  include  the  program  listings  for 
the  thermal  model,  the  sensitivity  program,  the  geometrical  preprocessing 
programs,  SCALC,  and  the  SRVC  absorption  model.  Also  included  in  the 
appendixes  are  the  geometrical  matrices  for  15  abstract  canopies,  the 
sensitivity  results,  and  supporting  validation  data. 


PART  II:  NEW  MODEL  STRUCTURE 


20.  This  part  summarizes  the  updated  formulation  and  solution 
approach  to  the  basic  thermal  canopy  model  developed  under  previous 
efforts.  The  individual  expressions  for  the  component  energy  budget 
processes  are  summarized  and  explicit  expressions  for  the  elements  of 
the  Jacobian  matrix  are  given.  The  geometrical  factorization  of  the 
energy  budget  equation,  particularly  for  the  longwave  flux  transfers, 
is  derived,  and  the  sequence  of  computer  programs  required  to  develop 
a  thermal  simulation  is  described. 

Energy-Balance  Framework 

21.  The  model  is  a  plane-parallel  abstraction  of  a  vegetation  canopy 
divided  into  three  horizontal  layers.  Two  additional  source  layers  are 
given  by  the  atmosphere  above  the  canopy  and  by  the  underlying  ground  or 
understory  layer.  An  energy-balance  framework,  assuming  steady-state 
conditions,  is  formulated  for  each  of  the  three  vegetation  layers  (sinks) 
as  a  function  of  the  five  source  layers.  For  this  and  subsequent  sections 
Figure  1  may  prove  useful  for  conceptualizing  the  various  energy  flows. 

The  sink  or  vegetation  layers  are  represented  by  i  =  1,2,3;  j=l,2,3,4,5 

represents,  respectively,  the  atmosphere,  the  three  vegetation  layers, 
and  the  ground  source  layers  of  energy  flux.  The  combination  of  the 

i , j  indices,  thus  represents  a  flow  of  energy  from  source  layer  j  to 
sink  layer  i  . 

22.  The  vector  expression  for  the  energy-balance  equations  was 
given  in  the  Part  I,  Equations  1  and  2  as: 


F  =  4  oo  B(X)TS  +  B(X)  +  A  +  H(X)  +  LE(X) 

23.  The  vector  equation  may  be  expanded  in  long  form  and  the  explicit 
dependence  on  parameters  or  input  variables  indicated  by 

%  a  jo  [B(Ta)Su  +  B(X1)S12  +  +  B(X3)S14  +  B(Tg)S15]  (8) 

+  Aj  -  aB(X1)  +  H(X1  ;WS,Ta)  +  LE(X1;WS,Ta ,R^ ,RH)  =  0 

H  a2°  [B<ys21  +  B(X1)S22  +  B(X2)S23  +  B(X3)S24  +  B{Tg)S25]  (9) 

+  A2  -  B(X2)  +  H(X2;WS,Ta)  +  LE(XZ;WS,T  .Rj.RH)  =  0 

4  [b(T,)S31  ♦  B(Xj)S32  ♦  B(X2)S33  ♦  B(X3)S34  *  B(Tg)S35]  (10) 

+  A3  -  aB(X3)  +  H(X3;WS,Ta)  +  LE(X3 ;WS ,Ta ,R] ,RH)  =  0 

where  the  explicit  formulation  for  each  energy  budget  component  used  in 


the  model  is  given  by 
Longwave: 


Bfx,)  =  E, (X,  +  273) 

»<V  *  ‘a<Ta  +  273>' 

B(Tg)  ■  Eg(T9  ♦  273)' 


Sensible  Heat:  H(X.WS,T)  =  (X,-TJ  -0.698(20.4  ♦  0. 2WSU-  ) ( 14) 

la  la 

Evapotranspi rati  on : 

LE(X - ;WS ,T  ,R, ,RH)  =  -697.75(-0.566  X.  +  597.3) 

l  a  I  1 

„  (5.234  e0,056715’Xi  -rh  5.234  e  0-056715-Ta)  10-6  /  ^ 

*  TTT9  V  *  ^  J 

R1  +  1/60  (0.04  +  1.27  WS  l/£) 


Shortwave  absorption:  A..  =  ABS(i)*SW 


where 


e  -  =  1  -  0.261  e'7*77  *  10  (17) 

ali 

ABS(i)  =  shortwave  absorption  coefficients  calculated  by  an  optical 
absorption  model  which  uses  a  Monte  Carlo  Technique  to 
include  multiple  scattering  effects  (see  Program  SRVC  in 
Appendix  A) 


Explicit  Evaluation  of  the  Jacobian 


24.  As  indicated  in  Part  I,  the  use  of  the  iterative  Newton- 
Raphson  technique  for  solving  the  nonlinear  thermal  equations  involves 
repeated  evaluation  of  the  expression 

6X  =  (jTJ)_1  JT  t-F(X0)]  (18) 


where : 

J  =  system  Jacobian  = 


(19) 


The  Newton-Raphson  method  is  employed  because,  in  this  case,  there  are 
relatively  simple  closed-form  expressions  for  the  elements  of  £  ,  and 
the  Jacobian  matrix  can  explicitly  be  evaluated.  Specifically, 


Ox 

3F2 

3Xi 

3Xt 

3X1 

!!i 

9F2 

sF3 

3X2 

3X2 

3X2 

9F1 

3F2 

9F3 

3X3 

3XS 

3X3 

The  i , j  component  of  J  is  easily  derived  as 


Jij  ■  2  «i  sj  sij  »<xj  *  273)3  *  5io  <4ej  "O'j  +  273)3 

(20.4  +  0.2WS0'97)  +  (697 .75) (0. 566) (5.234)( 10~6) (e0,056715 

R1  +  1/60  (0.04  +  1.27  WS 

-(697. 75 ) ( -0 . 566  X.  +  597. 3)(5. 234)(0. 056715)  10"6  e 
J  - 

R1  +  1/60  (0.04  +  1.27  WS  u'3) 
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+  0.698  T 

a 

XJ  -RH  e0-056715  Ta)  . 
-0. 5  j 

0.056715  Xi 

J 

(21) 


/ 


where : 

<5..  =  Dirac  delta  function 

25.  Program  TMODEL,  which  implements  the  equations,  is  given  in 
Appendix  A.  Subroutine  FEVAL  evaluates  the  function  and  the  Jacobian 
derivatives  and  calls  upon  Subroutine  BFUNC  which  calculates  the  long¬ 
wave  energy  component  and  derivative;  Subroutine  QFUNC  calculates  the 
sensible  heat  component  and  derivative.  It  should  also  be  noted  that  two 
different  expressions  for  the  convection  coefficient  arise,  depending  upon 
the  ambient  wind  speed.  Subroutine  RFUNC  calculates  the  evapotranspi ration. 

Geometrical  Factorization 

26.  A  significant  simplification  of  the  thermal  model  employed  in  this 
study  was  the  factorization  of  the  geometric-dependent  terms  from  the  energy- 
related  terms  for  the  longwave  flux  transfer  processes.  This  factorization 
is  made  possible  essentially  because  of  the  lack  of  multiple  scattering  in 

the  thermal  regime  between  canopy  components  whose  emissivities  (absorptivities) 
are  assumed  nearly  unity  and  by  the  fact  that  the  thermal  properties  on  both 
sides  of  a  canopy  component  are  assumed  equal.  The  significance  of  the 
factorization  is  not  so  much  in  the  increased  efficiency  in  model  calculation 
as  it  is  in  permitting  the  possibility  of  precalculating  these  geometrical 
matrices,  S  ,  for  a  wide  variety  of  plant  canopies.  These  precalculated 
matrices  may  then  be  convolved  with  the  appropriate  driving  variables  as 
required.  Program  SCALC  (Appendix  A)  performs  the  actual  calculations  for 
given  input  of  geometric  measurements. 

27.  The  required  input  data  for  a  three-layer  canopy  include 

f^  =  leaf  slope  distribution  for  layer  i= 1 ,2,3  and  angle 

G.  =5,15, . . .  ,85 
K  16 


* 


=  leaf  area  index  LAI  ,  for  layer  i 


Appendix  B  presents  the  matrices  calculated  for  five  different 
theoretical  canopies  at  three  different  LAI  densities  =1,4,  and  7. 

28.  The  five  theoretical  canopies  are  approximated  by  Verhoef 
and  Bunnik  (1975)  as 


Planopnile: 

Erectophile: 

Plagiophile: 

Extremophile: 

Uniform: 


f..  =  —  (1  +  cos  2  0.  ) 

ik  it  k' 

f...  =  -  (1  -  cos  2  0.  ) 

ik  ir  '  k' 

fik  =  f  (1  -  COS  4  6k) 

fik  =  I  (1  *  cos  4  0k) 


where  is  the  leaf  slope  angle. 

The  elements  of  the  S  matrix,  itself,  are  given  by 

9 
l 

k=l 


where: 


S..  =  E  f.,  C... 

ij  ik  i jk 


~-n 


2tt 


[ a •  r |  C0NT.jr  d«t>f  dOr 


(22) 


(23) 


a  is  the  orientation  of  the  leaf  at  angle  0k;  and  r  is 
the  direction  of  the  energy  flux  described  by  0  ,  <j> 


(i.e.,  r  =  (sin  0r  cos  $  ,  sin  sin  4>r,  cos  ©r)  (24) 


The  elements  of  CONT.^  represent  the  weighting  coefficients  which  give 
the  flux  contributions  from  a  source  layer,  j=l,2,3,4,5  ,  to  a  sink  vege¬ 
tation  canopy  layer,  i=l,2,3  ,  from  a  particular  source  direction  ©r,  $r  . 


/ 
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These  elements  for  an  arbitrary  direction,  r  ,  are  given  in  Table  1. 
PQ(i,r)  is  the  probability  of  a  gap  in  transversing  layer  i  at  direction 
r  .  It  may  approximated  by 

P0(1,r)  8  P0((,er)  •  e-"'1*  9<’-0r)  sec  (25) 

where  g ( i »0r )  is  the  mean  canopy  layer  projection  in  direction  . 

Mean  canopy  projection  is  given  by 


9(i.Qr)  = 


k(V\>  ftk  d8k 


where: 


k(Vk> 


=  2/tt  cos  0^  cos  ©r  ,  0^  i  ir/2  - 


=  4/tt  cos  0k  cos  Or  (<j>k  -  tt/2  -  tan<f>')  ,  Gk  >_  tt/2  -  ©r 


4>k  a  cos-1  (-cotQk  cot0r) 

Program  SCALC  also  calculates  the  view  factor  matrix  for  the  canopy. 

This  matrix  is  used  to  determine  the  thermal  flux  contribution  from  each 
vegetation  layer  and  the  ground  layer  which  is  intercepted  by  a  sensor 
viewing  the  canopy  at  a  particular  zenith  angle.  It  is  given  by 

W(i,0r)  =  VF(i,0r)  =  [VF(l.r)  VF(2,r)  VF(3,r)  VF(4,r)]T 


VF(l,0r)  =  1  -  PQ( 1 »r) 

VF(2,0r)  =  P0(l,r)[l  -  P0(2,r)l 
VF(3,0r)  =  PQ(l,r)  P0(2,r)[l  -  PQ(3,r )] 
VF(4,er)  =  P0(l,r)  P0(2,r)  PQ(3,r) 


c 


Sequence  of  Required  Computer  Runs 

29.  Appendix  A  contains  a  listing  of  all  the  computer  programs 
utilized  in  this  study.  Three  of  these  programs  are  directly  concerned 
with  thermal  modeling  or  preprocessing  steps  that  must  be  initiated  before 
the  thermal  calculations  may  be  made.  In  addition,  program  SENSIT  has 
been  included.  This  program  performs  the  systematic  and  repetitive  cal¬ 
culations  necessary  to  complete  the  sensitivity  calculations  of  many  of 
the  thermal  model  parameters. 

30.  The  basic  thermal  model  is  program  TMODEL.  This  program  assumes 
that  the  geometrical  characterization  of  the  canopy  has  been  performed  and 
the  appropriate  !B  matrix,  shortwave  absorption  vector,  and  view  factor 
matrix  have  been  calculated.  The  model  then  performs  similar  calculations 
at  discrete  time  intervals,  given  the  specification  of  the  appropriate 
parameter  (emission  and  absorption  characteristics  of  the  canopy  elements 
and  the  ground,  canopy  stomatal  resistance  to  water  vapor  diffusion). 
Furthermore,  the  input  information  must  be  provided  at  the  discrete  time 
intervals  simulated.  These  data  consist  of  the  air  temperature,  the  ground 
temperature,  the  wind  speed,  the  relative  humidity,  and  the  shortwave  flux. 
The  basic  philosophy  of  TMODEL  is  that  for  a  given  type  or  types  of  vege¬ 
tation  canopies,  one  would  want  to  simulate  a  multitude  of  scenarios  for 
their  thermal  behavior  based  on  either  ambient  meteorological  conditions 

or  modifications  to  the  thermal  properties  of  the  canopy  or  understory. 
Thus,  it  is  usually  required  to  calculate  the  geometrical  characteristics 
of  the  canopy  type  only  once  and  then  perform  multiple  simulations  of  the 
canopy  with  TMODEL. 

31.  The  calculations  of  the  appropriate  geometrical  flux  transfer 
matrices  are  done  by  Program  SCALC  and  Program  SRVC  for  absorption.  For 

both  of  these  programs,  detailed  canopy  geometry  information  is  required. 
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This  includes  the  leaf  area  index  for  each  layer,  and  the  leaf  slope 
distribution  by  layer.  In  addition,  to  calculate  the  shortwave  absorption 
coefficients,  average  optical  properties  of  the  canopy  elements  are  required. 
The  SRVC  absorption  model  is  further  described  in  a  report  by  Kimes,  Smith, 
and  Ranson  (1979). 

32.  The  complete  set  of  geometrical  matrices  have  been  calculated 
for  the  lodgepole  pine  canopy  in  Leadville,  Colorado,  studies  under  earlier 
WES  sponsorship,  the  Oouglas-fir  canopy  from  the  Cedar  River  Watershed, 
near  Seattle,  Washington,  and  the  oak-hickory  deciduous  community  at  the 
Walker  Branch  Watershed  at  Oak  Ridge  National  Laboratory  in  Tennessee.  In 
addition,  the  geometrical  characterization  has  been  performed  for  15 
abstract  canopies  of  varying  densities  and  geometrices.  These  data  are 
given  in  Appendix  B. 

33.  In  summary,  given  a  specific  canopy  to  be  studied  and  for  which 
detailed  geometrical  measurements  have  been  obtained,  Program  SCALC  and 
the  SRVC  absorption  model  are  first  used  in  a  preprocessing  manner  to 
calculate  the  appropriate  flux-transfer  matrices.  The  data  generated  from 
these  runs  are  then  used  in  Program  TMODEL.  If  there  is  no  specific  geo¬ 
metrical  measurement  available  for  canopies  of  interest,  then  one  of  the 
15  theoretical  canopies  in  Appendix  B  may  be  appropriate. 

34.  An  example  of  a  complete  analysis  for  the  validation  experiments 
is  given  in  Part  IV. 
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PART  III:  SENSITIVITY  ANALYSIS 


* 


35.  The  basic  analytic  model  described  in  this  report  may  be 
indicated  by  the  form: 


F(X,P,U)  =  0 

To  simplify  notation,  U  will  be  considered  to  be  an  additional  set  of 
parameters  augmenting  the  P  vector,  and  it  will  be  written  that: 

F(X,P)  =  0 


Further,  the  solution  to  the  system  of  equations  for  a  specific  parameter 

P„  will  be  indicated  as  X{P  ). 

-o  —  -o 

36.  Sensitivity  analysis  consists  of  determining  the  change  in  the 
solution  to  the  model  for  a  small  change  or  perturbation  in  model  parameters, 
i.e.,  XC^  +  AP). 

37.  The  sensitivity  function  S  is  defined  (Tomovic  1963)  as: 


1  im 
A  fM) 


+  AP)  -  X^) 
AP 


The  sensitivity  function  may  be  evaluated  analytically  by  differentiation 
of  the  system  equations  with  respect  to  the  parameters  under  consideration, 
yielding  the  following  sensitivity  equation: 


3F  3£  +  a£ 
al  3£  8£ 


or 


3F  +  3f 
3l  3P 


0 
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(27) 
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Alternatively,  computer  simulations  may  be  employed  in  which  the  parameters 
are  systematically  and  separately  perturbed  from  nominal  values  and  new 
canopy  temperatures  are  determined. 

38.  As  indicated  in  Part  I,  this  latter  approach  was  employed  for 
this  study.  Program  SENS  IT  was  written  to  facilitate  the  calculation 
(Appendix  A). 

39.  Program  SENSIT  requires  environmental  data  and  temperatures  for 
each  layer  to  initialize  the  analysis.  In  addition,  geometrical  factor 
matrices  describing  a  particular  canopy  are  required.  The  environmental 
data  used  was  collected  by  WES  personnel  at  Zweibriicken  Air  Force  Base  in 
West  Germany  on  4  October  1979.  Data  was  selected  at  0600  hours  and  1100 
hours  to  provide  for  nighttime  (predawn)  and  daytime  analysis.  Initial 
state  temperatures  for  each  layer  were  determined  from  simulation  results. 
The  sensitivity  analysis  was  performed  for  both  the  Douglas-fir  and  oak- 
hickory  canopies  resulting  in  a  total  of  four  analyses.  Table  2  lists  the 
initial  environmental  parameters  and  initial  temperatures  for  each  sensi¬ 
tivity  run.  Graphical  results  of  parameter  changes  versus  predicted  temp¬ 
eratures  are  found  in  Appendix  C. 

40.  The  daytime  sensitivity  analysis  showed  that  the  predicted 
canopy  temperatures  were  most  sensitive  to  the  air  temperature  within 
the  canopy.  A  10  percent  change  in  canopy  air  temperature  resulted  in 
nearly  a  10  percent  change  in  all  layers  for  both  types  of  canopies. 
Decreasing  longwave  absorption  coefficients  by  10  percent  resulted  in 
less  than  a  0.5°C  change  in  predicted  temperatures  and  showed  a  layer 
by  layer  dependence  for  both  canopies  and  time  periods.  Predicted 
canopy  temperatures  showed  minimal  sensitivity  to  changes  in  air  and 
ground  temperatures  as  input  to  the  model.  Temperature  predictions  were 

nearly  equally  sensitive  to  the  shortwave  absorption  in  all  three  layers 
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for  both  canopies.  Changing  the  canopy  emlssivity  in  the  top  layer  for 
both  canopies  had  little  effect  on  predicted  temperature  for  layer  1, 
but  slightly  increased  sensitivity  was  noted  for  the  two  lower  canopy 
layers.  Decreasing  ground  emlssivity  from  1.0  to  0.9  Increased  predicted 
temperatures  by  less  than  0.5°C.  Changing  relative  humidity  showed  little 
effect  on  canopy  temperatures  with  the  daytime  oak-hickory  analysis 
exhibiting  the  greatest  sensitivity.  A  linear  relationship  was  noted 
between  predicted  canopy  temperatures  and  the  parameters  discussed  above. 
Only  stomatal  resistance  and  wind  speed  analyses  showed  nonlinear  trends. 
Sensitivity  plots  of  stomatal  resistance  for  Douglas-fir  and  oak-hickory 
are  shown  in  Figure  2.  In  both  cases  the  plots  are  nonlinear  above 
Rj  values  of  0.08  min/cm.  Other  analyses  not  reported  here  showed  a 
linear  relationship  for  Rj  greater  than  0.08  min/cm  to  about  1.5  min/cm. 
Figure  3  shows  plots  of  wind  speed  versus  predicted  temperature  for  Doulgas 
fir  daytime  and  nighttime  analyses.  The  daytime  plot  shows  an  increase 
In  temperature  with  decreasing  wind  speed;  but  at  night,  temperatures 
decrease  slightly  with  decreasing  wind  speed. 
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PART  IV:  MODEL  VALIDATION  EXPERIMENTS 


c. 


41.  As  discussed  earlier,  the  objective  of  the  field  experiments 
was  to  provide  data  sets  from  diverse  targets  and  environmental  con¬ 
ditions  for  validation  of  the  Colorado  State  University  (CSU)  thermal 
canopy  model.  Two  existing. research  sites  were  located  through  the 
efforts  of  WES  personnel  that  proved  to  be  ideal  for  the  experiments. 

The  Cedar  River  site  was  located  in  a  Douglas-fir  forest  near  Seattle, 
Washington.  A  second  research  site,  the  Walker  Branch  Watershed,  was 
typical  of  an  Appalachian  deciduous  forest  and  was  located  near  Oak 
Ridge,  Tennessee.  Both  research  sites  were  being  used  for  ongoing 
research  in  forest  meteorology  and  possessed  extensive  instrumentation 
and  computerized  data  acquisition  support.  The  principal  scientist 
responsible  for  the  development  of  the  Cedar  River  site  was  Dr.  Leo  J. 
Fritschen  of  the  College  of  Forest  Resources,  University  of  Washington, 
while  Dr.  Boyd  A.  Hutchison  of  the  Atmospheric  Turbulence  and  Diffusion 
Laboratory  (ATDL) ,  National  Oceanic  and  Atmospheric  Administration 
(NOAA)  was  responsible  for  the  Walker  Branch  site.  Further  descriptions 
of  these  sites  are  given  below. 

Experimental  Design 

42.  The  model  validation  experiments  were  designed  by  CSU  and  WES 
personnel  with  cooperation  from  Drs.  Fritschen  and  Hutchison.  The  goal 
was  to  provide  appropriate  input  and  validation  data  for  the  CSU  canopy 
models.  Input  data  included  optical,  thermal,  and  environmental  para¬ 
meters  for  two  consecutive  24-hour  periods  of  the  targets.  Validation 


data  consisted  of  foliage  temperatures.  In  addition,  thermal  scanner 
imagery  was  to  be  obtained  by  local  National  Guard  units  at  specified 
times  throughout  the  measurement  periods.  Characterization  of  the 
foliage  angle  distributions  of  the  canopies  was  also  required.  Input 
data  requirements  and  methods  are  discussed  in  a  later  section. 

43.  WES  personnel  were  responsible  for  overall  mission  coordination, 
thermal  radiometric  measurements  of  ground  and  canopy,  air  temperature 
measurements  in  the  lower  1.5  m  of  the  canopies,  and  arranging  for 
National  Guard  thermal  scanner  overflights  of  the  experimental  sites. 

CSU  personnel  comnunicated  requirements  for  micrometeorological  data 
to  Drs.  Fritschen  and  Hutchison,  obtained  foliage  geometry  data  from 
the  sites,  and  performed  necessary  optical  measurements  required  to 
run  the  canopy  models.  Groups  headed  by  Drs.  Fritschen  and  Hutchison 
provided  site  access,  operated  and  maintained  the  data  acquisition 
systems,  and  provided  assistance  for  interpreting  the  micrometeorological 
data.  In  addition.  Dr.  L.  W.  Gay  of  the  School  of  Renewable  Natural 
Resources,  Arizona  University  at  Tucson  participated  in  the  Cedar  River 
Douglas-fir  experiment  to  test  the  use  of  direct  beam  depletion  measure¬ 
ments  for  determining  forest  biomass. 

Site  Descriptions 

44.  Two  established  research  sites  were  available  for  this  study. 

A  site  near  Seattle,  Washington,  developed  and  maintained  by  Dr.  Leo 
Fritschen  of  the  University  of  Washington,  provided  data  for  a  stand 
of  mature  Douglas-fir.  Dr.  Boyd  Hutchison  of  ATDL/NOAA  made  available 
an  oak-hickory  site  near  Oak  Ridge,  Tennessee  and  provided  necessary 
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environmental  data.  A  detailed  description  of  these  sites  is  provided 
below. 


Cedar  River,  Washington 

45.  The  Cedar  River,  Washington,  study  site  is  located  on  the  A.  E. 
Thompson  Research  Center  at  the  western  end  of  the  Cedar  River  Watershed. 
The  site  lies  in  the  Puget  Sound  Basin  at  the  western  foot  of  the  Cascade 
Mountains  55  km  southeast  of  Seattle,  Washington,  at  47°23'N  and  121° 

56 ' W .  The  elevation  is  approximately  215  m  above  mean  sea  level. 

46.  The  area  was  logged  prior  to  1924  and  subsequent  fires  resulted 
in  a  mosaic  of  different  aged  stands  (Jensen,  1976).  The  most  common 
community  on  the  site  is  Douglas-fir  (Pseudotsuqa  menziesii  (Mirb.) 
Franco).  This  naturally  regenerated  stand  was  approximately  41  years 
old  with  an  average  tree  spacing  of  5.8m.  There  were  572  trees  per 
hectare  consisting  mainly  of  Douglas-fir,  a  few  hemlock,  and  maple 
(Figure  4).  Ground  cover  consisted  of  fern,  salal,  huckleberry  mosses, 
and  litter  (Figure  5).  Bare  soil  areas  were  minimal  and  occurred  only 

on  roads  and  other  localized  disturbed  areas.  Soil  at  the  site  con¬ 
sisted  of  Barneston  gravelly,  1 oamy  sand  originating  from  glacial 
outwash. 

47.  The  specific  study  site  was  located  at  a  micrometeorological 
observatory  maintained  and  operated  by  the  University  of  Washington. 
Average  height  of  the  Douglas-fir  stand  was  about  28  m  with  an  average 
LAI  of  approximately  7.8.  Located  at  this  site  was  a  28-m-tall  Douglas- 
fir  tree  contained  in  a  lysimeter  (Fritschen,  Cox,  and  Kinerson,  1973). 
The  site  adjacent  to  this  tree  was  instrumented  to  provide  data  for 
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evapotranspi ration  studies.  These  data  included  wet  and  dry  bulb 
temperatures,  soil  temperatures,  global  shortwave  radiation,  pre¬ 
cipitation,  and  wind  speed  and  direction.  In  addition,  needle  surface 
temperatures  were  monitored  at  several  points  around  the  lysimeter  tree 
near  the  top  and  center  of  the  canopy.  These  data  were  recorded  at 
selected  time  intervals  by  a  computerized  data  acquisition  system. 

A  33-m  walk-up  tower  was  available  adjacent  to  the  lysimeter  tree  to 
provide  access  to  needle  temperature  sensors  and  other  measurement 
devices. 

Walker  Branch,  Tennessee 

48.  The  Walker  Branch  study  site  is  located  near  the  Walker  Branch 
Watershed  research  facility  on  the  U.S.  Department  of  Energy  Reservation 
near  Oak  Ridge,  Tennessee,  at  35°58'N  and  84°15'W.  An  intensive  forest 
meteorological  research  site  operated  by  the  ATDL  of  the  NOAA  was  made 
available  for  this  study.  This  research  area  is  situated  on  a  ridge 

top  about  70  m  above  the  valley  floor  at  an  elevation  of  335  m  above  mean 
sea  level. 

49.  The  area  is  representative  of  an  Appalachian  deciduous  forest 
(Hutchison,  1977).  The  species  composition  of  the  stand  is  dominated 
by  various  species  of  oak  and  hickory,  including  Quercus  alba,  Quercus 
prinus,  Quercus  vel uti na ,  Carya  glabra  and  Ca.rya  ovata.  Acer  rubra  (red 
maple),  Prunus  serotina  (black  cherry),  Liriodendron  tulipifera  (yellow 
popular)  are  less  frequently  found.  Common  understory  plants  include 
Oxydendron  arboreum  (sour  wood),  Cornus  florida  (flowering  dogwood) 

and  Cercis  canadensis  (eastern  redbud).  The  average  height  of  the 
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codominant  trees  forming  the  canopy  is  about  21.5  m  with  lower  limit 

of  the  live  crown  being  15  m  above  the  ground.  These  heights  vary 

greatly  due  to  the  uneven  age  of  the  stand  (Figure  6).  Basal  area  was 
2  -1 

approximately  26  m  ha  .  The  site  appeared  parklike  due  to  a  fire  that 
occurred  several  years  ago.  Understory  growth,  however,  is  abundant. 

The  ground  is  covered  by  an  accumulation  of  litter  (Figure  7)  with 
bare  soil  occurring  only  in  disturbed  areas.  In  addition,  fragmented , 
grey-colored  rock  covered  the  road  surfaces.  A  metal  track  was  in  place 
beneath  the  stand  to  provide  all-weather  access  for  research  vehicles. 

This  track  was  covered  with  litter  by  ATDL  personnel  during  the  field 
experiments . 

50.  The  site  is  extensively  instrumented  to  record  data  pertinent 
to  forest  meteorology  research  as  well  as  the  thermal  modeling  studies. 
Hutchison  (1977)  gives  a  detailed  description  of  the  research  facility. 

Modeling  Input  Data 

51.  The  data  collected  at  the  two  sites  included  foliage  and  back¬ 
ground  optical  parameters,  geometry  characterization  measurements ,  and 
environmental  measurements.  This  section  describes  the  data  required 

for  the  models  and  the  techniques  or  sources  used  to  acquire  it.  Listings 
of  the  data  values  are  included  in  Appendix  D. 

Foliage  geometry 

52.  The  structure  of  a  canopy  defined  by  the  foliage  inclination 
angles  and  LAI  is  important  for  characterizing  the  interactions  of 
radiation  with  the  canopy.  These  inputs  are  required  by  the  optical  SRVC 
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model  (Oliver  and  Smith  1974)  to  estimate  the  shortwave  absorption  of 
a  canopy  and  by  the  thermal  model  to  describe  longwave  energy  exchanges 
inside  and  outside  of  the  canopy. 

53.  The  procedure  for  determining  foliage  geometry  included  acquiring 
high-contrast  black-and-white  slide  photography  of  canopy  silhouettes. 

These  slides  serve  as  input  to  a  laser  diffractometer  which  characterizes 
the  frequency  of  occurrence  of  foliage  angles  in  terms  of  the  resulting 
diffraction  pattern.  The  diffraction  patterns  are  optically  sampled,  and 
the  results  are  analyzed  with  a  series  of  computer  programs.  See  Kimes, 
Smith,  and  Ranson  (1979)  for  a  discussion  of  the  theory  and  procedures. 

54.  The  walk-up  towers  at  both  sites  provided  an  excellent  platform 
for  acquiring  slides  of  the  canopies.  For  the  purposes  of  the  modeling, 
the  canopies  were  partitioned  into  three  layers  of  equal  height.  Photo¬ 
graphs  were  taken  for  each  layer  from  several  directions  from  the  tower. 
This  provided  a  larger  sample  size  and  minimized  effects  of  azimuthal 
asymmetry.  Ideally,  the  photographs  should  be  taken  with  a  white  back¬ 
drop  placed  behind  the  target  to  eliminate  background  trees  and  shadows. 
However,  this  was  impractical  for  the  canopies  under  study.  As  a  result, 
the  slides  were  manually  interpreted  to  delineate  branches  of  the  desired 
tree  in  the  photographs.  This  was  done  by  projecting  the  slide  on  white 
paper  and  tracing  the  appropriate  branches.  Earlier  work  by  Kimes, 

Smith,  and  Ranson  (1979)  showed  that  for  complex  canopies,  such  as 
conifers,  two  interpretations  are  required:  one  with  all  branches 
represented,  and  a  separate  tracing  including  only  branches  bearing 
foliage.  High-contrast  slides  of  these  tracings  were  used  as  input  to 
the  laser  diffractometer.  The  branch  and  foliage  measurements  were 


combined  later  to  provide  the  inclination  angle  distributions  for 
each  layer. 

55.  The  calculated  foliage  angle  distributions  for  a  Douglas-fir 
canopy  are  shown  in  Figure  8.  For  comparison  purposes,  distributions 
of  lodgepole  pine  (Pinus  contorta)  reported  by  Kimes,  Ranson,  Kirchner, 
and  Smith  (1978)  are  included.  Figure  9  shows  foliage  angle  distributions 
for  oak-hickory.  These  data  were  derived  from  direct  measurements  pro¬ 
vided  by  Dr.  Hutchison.  Laser  diffraction  results  for  oak-hickory  were 
unavailable  due  to  equipment  problems.  For  comparison  a  one-layer  dis¬ 
tribution  for  Russian  olive  (Elaeaqnus  anqusti fol ia)  reported  by  Kimes, 
Smith,  and  Ranson  (1979)  is  included. 

Leaf  area  index 

56.  LAI  is  defined  as  the  total  one-sided  leaf  area  occupying  the 
horizontally  projected  area  of  the  canopy.  For  example,  an  LAI  of  5 
indicates  that  five  layers  of  leaves  could  be  overlayed  to  completely 
fill  an  area  equal  to  the  canopy  projection  on  the  ground.  LAl's  for 
this  study  were  determined  from  data  provided  by  Drs.  Fritschen  and 
Hutchison.  LAl's  for  the  Douglas-fir  canopy  were  derived  from  measure¬ 
ments  reported  by  Kinerson  and  Fritschen  (1971).  In  this  report,  graphs 
of  canopy  height  .z(m)  versus  surface  area  density  F(_z)  (m  m  )  for 
nine  sample  plots  are  given.  Integrating  £(,z)  over  height  gives  the 
needle  surface  area  index  NSAI  for  a  particular  height  increment  dz. 
Data  points  were  taken  from  the  graphs  and  averaged  for  given  heights  to 
produce  a  single  average  surface  density  curve.  This  curve  was  parti¬ 
tioned  into  three  layers  of  equal  height  and  layer  NSAI's  determined 
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by  Simpsons  Rule  (Figure  10).  For  our  modeling  purposes,  LAI  values 
were  determined  by  dividing  NSAI  for  each  layer  by  two. 


57.  LAI  for  the  oak-hickory  canopy  was  determined  from  data 
provided  by  ATDL.  These  data  consisted  of  a  graph  of  cumulative  LAI 
versus  height  and  graph  of  LAI  at  given  heights  through  the  canopy. 

A  smoothed  version  of  the  latter  is  presented  as  Figure  11. 

Canopy  density  parameter 

58.  This  parameter  ranges  from  0  to  1  and  describes  the  spatial 
dispersion  of  foliage  elements  within  a  canopy.  As  values  approach  1, 
gaps  in  the  canopy  are  less  frequent  since  the  foliage  is  more  regularly 
dispersed.  This  parameter  is  used  in  the  equation  to  determine  the 
probability  of  gaps  occurring  in  a  canopy  layer.  A  value  of  0.1  was 
chosen  for  all  model  runs.  For  a  detailed  discussion  of  spatial  dis¬ 
persion  of  canopies  see  deWit  (1965). 

Canopy  optical  paran-  ^ ers 

59.  The  shortwave  transmission  and  reflectance  of  foliage  elements 
are  required  as  inputs  for  estimating  average  absorption  coefficients 

as  discussed  below.  Canopy  element  transmission  values  were  measured 
at  the  study  sites,  but  reflectance  values  were  derived  from  the 
published  literature. 

60.  The  procedure  for  determining  transmission  consisted  of  placing 
a  needle  or  leaf  over  a  narrow  slit  on  a  flat  plate  attached  to  a  photo¬ 
diode  and  recording  a  reading  of  the  amount  of  light  passing  through  the 
sample.  Measurements  were  made  in  four  wavelength  bands--at  4.8um, 
0.55pm,  0.68pm,  and  0.80pm.  The  transmission  measurements  were  then 
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ratioed  to  the  Incoming  spectral  irradiance  measured  from  a  BaSo^ 
standard  reflectance  panel.  The  measurement  procedure  was  repeated 
for  several  foliage  samples  and  the  results  averaged.  Natural  illumi¬ 
nation  was  used  for  the  Douglas-fir  needles;  however,  because  of  rapidly 
changing  irradiance  conditions  at  the  Walker  Branch  site,  a  bank  of 
fluorescent  tubes  was  used  as  the  irradiance  source.  The  transmission 
measurements  were  integrated  over  wavelength  to  estimate  the  average 
shortwave  transmittance  from  0.48  to  0.80pm.  This  wavelength  interval 
was  assumed  adequate. 

61.  Shortwave  reflectance  values  for  Douglas-fir  were  obtained  from 
data  presented  by  Jarvis,  James,  and  Landsberg  (1976).  Curves  for  old 
and  new  Douglas-fir  needles  were  digitized  and  averaged.  The  resulting 
curve  was  then  integrated  over  the  wavelength  interval  from  0.45pm  to 
1.2pm  to  obtain  the  average  shortwave  reflectance  coefficient.  The  oak- 
hickory  canopy  element  reflectance  was  determined  from  data  presented 

by  Colwell  (1969).  Data  for  maple,  oak,  and  yellow  poplar  were  averaged 
and  integrated  over  the  wavelength  interval  0.45pm  to  1.2pm. 

62.  In  addition  to  foliage  transmission  and  reflectance  estimates, 
an  average  background  reflectance  was  determined  at  both  sites.  Measure¬ 
ments  were  made  of  various  surface  covers  such  as  litter,  bare  soil,  and 
ground  cover  vegetation.  The  results  were  weighted  according  to  visual 
estimates  of  occurrence  and  then  averaged  and  integrated. 

Shortwave  absorption  coefficients 

63.  The  absorption  of  global  shortwave  radiation  by  canopy  layers 

is  an  important  component  in  the  daytime  energy  budget.  It  is,  however, 

difficult  to  directly  measure  and  must  be  estimated  with  models.  These 
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coefficients  were  approximated  with  the  SRVC  model  modified  for  absorp¬ 
tion  (Kimes,  Smith,  and  Berry,  1980).  The  procedure  involved  running 
model  simulations  with  appropriate  canopy  layer  geometry,  LAI  ,  and 
optical  parameters  for  an  average  zenith  sun  angle  of  45°.  The  resulting 
absorption  values  represent  the  proportion  of  shortwave  absorption  in  each 
canopy  layer.  Since  the  thermal  model  requires  absorption  per  unit  leaf 
area,  the  simulated  absorption  coefficients  were  divided  by  the  one-sided 
leaf  area  in  a  given  layer. 

Stomatal  resistance 

64.  The  resistance  of  the  leaf  to  water  vapor  diffusion  depends  on 
many  environmental  factors.  Leaf  stomates  open  and  close  in  response 
to  microclimatic  and  soil  conditions  and  regulate  the  cooling  of  the 
plant  through  evapotranspi ration .  Thus,  stomatal  resistance  is  important 
when  considering  energy  budget  analysis  of  plants.  This  parameter  is 
difficult  to  measure,  so  for  modeling  purposes  average  values  were  used 
as  constants.  The  value  for  Douglas-fir  was  set  at  0.66  min/cm  as  an 
average  value  for  coniferous  forest  (Kimes,  Smith,  and  Ranson,  1979). 
Stomatal  resistances  were  determined  from  data  provided  by  Hutchison*. 
These  data  ranged  from  0.04  to  0.07  min/cm  for  sun  leaves.  The  upper 
value  was  selected  for  use  in  all  deciduous  canopy  simulations. 


*  Personal  communication;  B.  A.  Hutchison,  Atmospheric  Turbulence  and 
Diffusion  Laboratory,  National  Oceanic  and  Atmospheric  Administration, 
Oak  Ridge,  Tennessee,  1979. 
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Emissivity  and  absorptivity 

65.  The  ability  of  a  canopy  element  to  emit  and  absorb  longwave 
radiation  is  expressed  by  the  emissivity  and  absorptivity  coefficients 
specified  for  each  component  in  the  canopy  layers  and  for  the  ground 
layer.  Available  literature  values  or  direct  measurements  could, 
consequently,  be  substituted.  For  all  of  the  analyses  reported  here, 
the  emissivity  and  absorptivity  a..  are  set  equal  to  1.0  for 

each  of  the  three  canopy  layers.  Emissivity  of  the  ground  e^  was  also 
set  at  1.0.  Emissivity  of  the  air  e5  was  calculated  as  a  function 

a 

of  air  temperature  by  the  following  function  (Hudson,  1969): 


e,  =  1.0  -  0.0261  e(-°-000777  Ta  ) 

a 


Canopy  Temperature  Measurements 

66.  Since  the  purpose  of  the  experiments  was  to  collect  data  sets 
for  validation  of  the  thermal  model,  actual  canopy  foliage  temperature 
measurements  were  required.  The  experiments  were  designed  to  provide 
measured  canopy  temperatures,  as  well  as  thermal  scanner  images  of 

the  sites. 

67.  The  experimental  setup  at  the  Cedar  River  site  included  tem¬ 
perature  measurements  for  a  number  of  individual  Douglas-fir  needles.  The 
temperature  sensors  were  located  around  the  lysimeter  tree  at  average 
heights  of  26  m  and  20  m.  The  measurements  at  a  given  height  were 
averaged  to  give  an  average  layer  measurement.  The  26-m  measurement 

was  assumed  to  represent  the  average  canopy  temperature  for  layer  1. 

The  20-m  measurement  approximated  layer  2,  although  its  location  was 
closer  to  the  boundary  between  layer  1  and  layer  2.  These  layer  temp¬ 
eratures  are  plotted  along  with  air  temperature  against  time  in  Figure  12. 
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68.  No  individual  leaf  temperature  measurements  were  available  at 
the  Walker  Branch  site,  so  a  portable  thermal  radiometer*  was  used  to 
monitor  the  canopy  temperature  throughout  a  24-hour  period.  The  pro¬ 
cedure  was  to  position  the  instrument  upward  from  the  ground  at  the 
canopy  and  slowly  move  it  until  the  maximum  temperature  was  recorded. 

This  was  done  to  minimize  errors  due  to  the  presence  of  sky  or  clouds 
in  the  field  of  view.  Figure  13  shows  a  plot  of  the  canopy  temperature 
with  air  temperature  above  the  canopy  and  ground  temperature  against 
time. 

69.  In  addition  to  the  geometrical,  optical,  and  thermal  parameters 
discussed  above,  a  set  of  dynamic  variables  characterizing  the  micro¬ 
climate  of  the  target  is  required  to  drive  the  thermal  model.  These 
parameters  consist  of  air  temperature  above  the  canopy,  ground  surface 
temperature,  wind  speed  at  the  top  of  the  canopy,  relative  humidity,  and 
global  shortwave  radiation. 

70.  Air  temperature,  ground  temperature,  and  shortwave  radiation 
are  important  components  for  energy  exchange  into  and  within  the  system; 
whereas  wind  speed  and  relative  humidity  are  important  for  determining 
forced  convection  loss  and  evapotranspirative  cooling  of  plants, 
respectively. 

71.  Environmental  data  were  provided  from  the  automated  recording 
systems  at  the  two  sites.  Air  and  ground  temperatures  and  global 
shortwave  radiation  were  measured  directly.  Relative  humidity  was 
determined  from  wet  and  dry  bulb  temperatures.  All  measurements  were 

*  Barnes  Insta-Therm,  Barnes  Engineering  Corporation. 
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either  instantaneous  or  short  time  interval  averages.  Plots  of  the 
four  environmental  parameters  are  shown  in  Figure  14  for  Cedar  River 
and  Figure  15  for  Walker  Branch. 
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Model  Validation  Results 

72.  The  data  collected  for  the  coniferous  Doulgas-fir  and  deciduous 
oak-hickory  canopies  provided  a  good  means  of  testing  the  thermal  model 
under  these  diverse  conditions.  Three-layer  canopy  temperature  simu¬ 
lations  were  made  over  a  48-hour  period  with  both  data  sets  and  the 
results  were  compared  with  measured  temperatures . 

Douglas-fir  canopy 

73.  The  thermal  model  was  run  with  environmental  data  acquired  over 
the  48-hour  period  of  4-5  August  1979.  These  data  plus  the  required 
geometrical  factor  matrices  which  include  the  longwave  exchange  co¬ 
efficients,  the  sensor  view  angle  weighting  factors,  and  average  short¬ 
wave  absorption  coefficients  are  listed  in  Appendix  D.  The  emissivities 
and  thermal  absorption  coefficients  for  each  layer  were  set  to  1.0.  The 
total  canopy  resistance  to  water  vapor  diffusion  was  input  at  0.66  min/cm. 

74.  A  plot  of  the  simulated  three-layer  temperatures  with  measured 
air  temperature  is  shown  ir.  Figure  16.  The  layer  1  simulated  temperatures 
follow  the  trend  of  air  temperature,  but  fall  below  during  the  night  and 
are  higher  during  the  day.  The  layer  2  and  layer  3  predictions  are 
nearly  equal  to  air  temperature  throughout  the  48-hour  period.  Com¬ 
parisons  of  measured  and  predicted  needle  temperatures  for  layers  1  and 

2  are  presented  as  Figures  17  and  18,  respectively. 
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75.  The  layer  1  predicted  temperatures  vary  from  the  measured 
temperature  by  a  maximum  of  3°C.  These  deviations  were  observed  during 
the  daytime  hours  under  very  hazy  skies.  Nighttime  predictions  deviated 
from  measured  by  2°C  or  less  with  the  maximum  deviations  occurring  under 
conditions  of  fog.  This  leads  to  the  conclusion  that  the  thermal  model 
may  be  most  valid  for  days  with  primarily  direct  solar  radiation  and 
clear  nights  where  radiative  cooling  is  occurring. 

Oak-hickory  canopy 

76.  Environmental  data  acquired  at  the  Walker  Branch  site  for  the 
48-hour  period  from  18-19  August  1979  were  used  to  validate  the  thermal 
model  for  a  deciduous  oak-hickory  canopy.  Emissivities  and  thermal 
absorption  coefficients  for  the  three  canopy  layers  were  set  to  1.0. 
Canopy  resistance  to  water  vapor  diffusion  was  input  as  0.07  min/cm 
and  held  constant.  The  input  environmental  and  geometrical  factor  data 
for  this  canopy  simulation  are  presented  in  Appendix  D. 

77.  Figure  19  presents  the  three-layer  canopy  temperature  pre¬ 
dictions  along  with  measured  air  temperature.  Nighttime  simulations 
were  nearly  equal  to  air  temperature,  but  daytime  predictions  varied 
by  a  maximum  of  2°C  over  air  temperature. 

78.  Measured  temperatures  were  compared  to  predicted  results  for 
layer  2  and  are  shown  in  Figure  20.  The  agreement  between  model  and 
measured  temperatures  was  quite  good.  The  largest  deviation  (3°C) 
occurred  in  the  afternoon;  but  morning  and  nighttime  predictions  varied 
by  only  1°C  or  less. 
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Summary 

79.  The  results  of  the  model  validation  study  indicate  that  the 
thermal  canopy  temperature  model  provided  good  estimates  of  actual 
temperatures  for  nighttime  periods  to  within  2°C  for  both  canopies 
studied.  Daytime  simulations  generally  underestimated  measured  tem¬ 
peratures  for  Douglas-fir  and  overestimated  temperatures  for  oak-hickory. 
The  results  indicate  that  the  model  may  not  adequately  account  for 
energy  transfers  under  foggy  or  very  hazy  conditions. 


PART  V: 


RECOMMENDATIONS 


80.  Two  broad  directions  for  further  research  and  development  are 
suggested  in  the  paragraphs  below.  The  first  set  of  tasks  represent 
logical  extensions  or  improvements  to  the  thermal  model  utilized  in  this 
study.  Also,  a  not-quite-so-obvious  extension  to  the  calculation  of  the 
geometrically  dependent  flux  transfer  matrices  is 'outlined.  The  second 
thrust  recommended  for  further  development  is  concerned  with  parameter 
estimation  techniques  which  can  be  used  to  estimate  model  parameters, 
control  (or  input)  variables,  and  elements  of  the  state  vector  itself. 

Two  techniques  are  described.  The  first  technique  based  on  sensitivity 
functions  is  appropriate  for  the  steady-state  version  of  the  model.  The 
second  method,  based  on  the  Kalman  filter,  is  more  appropriate  for  dynamic 
representation  of  the  thermal  model . 

Model  Improvements 

81.  The  most  urgent  need  for  model  improvement  is  to  evaluate  differ¬ 
ent  theories  for  the  height  dependence  of  air  temperature  within  the  canopy 
and  of  the  vertical  profile  for  wind  speed.  It  is  particularly  appropriate 
to  examine  those  techniques  which  would  yield  these  temperature  and  wind 
profiles  from  a  few  limited  measurements.  The  structure  of  the  current 
thennal  model  can  easily  include  vertical  variations  in  the  two  parameters; 
they  are  held  constant  for  the  want  of  better  knowledge  and  for  simplicity. 

82.  The  utility  of  the  model  could  be  extended  if  a  ground  temperature 
module  was  included.  Particularly  for  this  extension  it  may  be  appropriate 

to  develop  a  time-dependent  version  of  the  model  to  include  heat  storage  effects. 
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83.  A  useful  exercise,  but  of  lesser  priority,  would  be  to  systematically 
examine  the  alternative  formulations  expressed  by  various  authors  for  dif¬ 
ferent  components  of  the  energy  budget  equation;  that  is,  evapotranspi ration, 
sensible  heat,  and  so  forth.  There  is  no  clear  rationale  for  selecting  one 
expression  over  another.  However,  the  separate  expressions  can  be  programmed 
and  sensitivity  analysis  performed  on  the  individual  expressions. 

84.  Finally,  further  analysis  of  the  structure  of  the  geometrical 
matrices  should  be  carried  out  to  determine  if  either  an  analytical 
decomposition  of  the  matrices  into  a  leaf  density  (leaf  area  index)  component 
and  leaf  slope  distributions  can  be  constructed.  If  an  analytical  decompo¬ 
sition  is  not  possible,  then  numerical  interpolation  techniques  should 

be  investigated. 

85.  As  an  example,  consider  the  expressions  for  the  view  factor 
matrix  VF(i  ,0)  where  0r  is  the  zenith  view  angle  and  i=l,2,3,4 
corresponds  to  contributions  from  the  three  vegetation  layers  and  the 
ground  surface: 

VF(i,or)  =  [VF(l,r)  VF(2,r)  VF(3,r)  VF(4,r)]T 
VF(l,0r)  =  1  -  PQ(l,r) 

VF(2,0r)  =  PQ(l,r)  (1  -  Po(2,r) 

VF(3,0r)  =  PQ(l,r)  Po(2,r)  (1  -  PQ(3,r)) 

VF(4,er)  =  PQ(l,r)  Po(2,r)  PQ(3,r) 

where: 

Po(1,0>  =  e'LAI  9(f'0)  sec  0 

LAI  (i)  =  the  mean  leaf  area  index  for  layer  i 

g(i,o)  =  the  mean  canopy  projection  of  vegetation  layer  i  in 

the  direction  0  ,  depending  only  on  the  leaf  slope  dis¬ 
tributions  for  layer  i 
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86.  A  direct  factorization  is  not  apparent.  However,  particularly 
for  large  LAI  a  Taylor  series  expansion  would  yield  a  more  tractable 
form.  Alternately,  LAI  could  be  varied  between  0  and  10  and  numerical 
tables  generated. 

Parameter  Estimation 

87.  Two  different  approaches  are  suggested  for  estimation  of  param¬ 
eters,  control  vector  inputs,  and/or  selected  components  of  the  unknown 
state  vector,  that  is  the  average  canopy  temperature  for  the  three  different 
layers.  One  approach  is  more  applicable  to  the  steady-state  conditions; 

the  second  approach  is  more  appropriate  for  the  time-dependent  version  of 
the  model.  In  each  case  it  is  assumed  that  selected  measurements  of  canopy 
temperatures  are  available  for  some  time  periods  and  that  some  of  the 
parameters  and  control  vector  components  are  also  known.  A  typical  scenario 
would  be  that  the  top  layer  canopy  temperature  is  measured  over  a  diurnal 
cycle  and  that  all  parameters  and  input  components  are  known  except  for 
the  ^  matrix,  the  longwave  flux  transfer  matrix.  It  is  then  desired  to 
estimate  the  S  matrix  which  depends  on  the  geometrical  properties  of  the 
canopy  and  evaluate  the  fit  on  a  second  diurnal  cycle.  Other  scenario 
examples  can  be  envisioned.  In  this  section,  general  development  of  the 
two-parameter  estimation  techniques  are  indicated. 

88.  First,  consider  the  steady-state  situation  where  the  model  is 
given  by  the  following  equation: 

F(X,P,U)  «  0 

where  the  symbols  have  the  same  meaning  as  given  earlier.  For  this 
situation  the  parameter  estimation  technique  of  nonlinear  systems  as 
described  by  Durando  and  Leondes  (1976)  is  recommended.  For  simplicitity 
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the  U  vector  is  appended  to  the  £  vector  and  the  equation  is  re¬ 
expressed  as: 

F(X,P)  =  0 


Further,  it  is  assumed  that  observation  variables  are  the  canopy  temperature 
variable,  x  ,  directly.  Given  a  known  measurement,  becomes 

a  function  of  only.  Assume  an  initial  estimate  of  F\P^  .  Then 
£(X^PJ  can  be  expanded  about  PQ  : 

£<*.p>  -  f<vV  ■  |.p  (p-V  +  ^  (28) 

—  --o 

For  the  steady-state  formulation  £(X,P)  =  0  ;  e  is  the  error  vector. 
Iteration  is  continued  until  convergence,  i.e., 


^n+I 


P  +  6P 

— n  — 


(29) 


89.  If  observations  are  available  for  more  than  one  time  interval, 
the  optimal  P^  is  chosen  which  minimizes  the  sum  of  eTe  over  all  time 
intervals.  More  general  formulations  of  this  approach,  including  the  use 
of  a  variable  increment  step  size,  are  given  in  the  paper  by  Durando  and 
Leondes. 

90.  The  second  technique  proposed  is  applicable  to  the  time-dependent 
formulation  of  the  thermal  model  given: 

9X 

M  -r  *  F(X,P,U,T)  (30) 


where: 

M  =  specific  heat  capacity  of  the  system 
T  =  time 
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The  general  approach  recommended  here  is  the  use  of  the  Kalman  filter 
after  first  linearizing  the  system.  Specifically, 


X  =  AX  +  BU  +  W 

Z  =  H  X  +  V 


(31) 

(32) 


where  X.  =  3X/3T  represents  the  dynamical  equations  of  the  system,  A 
and  B  are  expansion  matrices,  and  W  represents  the  modeling  error. 

91.  1  is  the  observation  vector,  which  now  permits  transformation 

on  the  state  vector  (canopy  temperature),  and  V  is  the  observation  noise. 

92.  Kalman  filtering  on  the  state  vector  or  on  the  augmented  state 
vector,  that  is,  after  appending  or  U  to  X.  »  1S  then  given  by  the 
standard  expressions  (Friedland  1972): 


where: 


X  =  X  +  K  (ZB  -  H  XJ 
n  n  n  n  n 


Xn  =  *n  .  Xn  . 
n  n-l  n-l 


K  =  P  HT  (H  P  HT  +  V  )_1 
n  n  n  '  n  n  n  ny 


^n  =  $n-l  ^n-1  n-l  +  Bn-1  Wn  Bn-1 


P  *  (I  -  K  H  )  P 
n  n  n  n 


(33) 

(34) 

(35) 

(36) 

(37) 


4>  is  the  transition  matrix  for  the  system,  n  represents  the  discrete 
time  interval,  and  X  describes  the  model  predictions. 

93.  An  additional  is  required  if  many  time  intervals  are 
available,  e.g.,  a  diurnal  cycle;  however,  the  final  estimates  are  in¬ 
sensitive  to  these  values. 
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Table  1 


Table  2 


Initial  environmental  and  initial  temperature  data  used 
for  sensitivity  analyses  for  the  Douglas-flr  and  oak-FTTckory 
canopies 


Enviro nmental  Data 


Time 

Jr 

qT 

WS 

swr2 

hours 

°_C 

cm/s 

RH 

w/nr 

0600 

10.6 

10.7 

136.0 

0.72 

0.0 

1100 

18.2 

19.0 

110.0 

0.84 

299.7 

Initial  Temperatures,  °C 
Time 

hours  Layer  1  Layer  2  Layer  3 


Douglas-fir 

0600 

9.0 

10.1 

10.1 

1100 

18.4 

18.2 

18.2 

Oak-hickory 

0600 

10.1 

10.5 

10.5 

1100 

18.8 

18.5 

18.2 

/ 


SINK 

VARIABLE 


SOURCE 

VARIABLE 


SKY  J=1 

1-1  VEGETATION  LAYER  ONE  J=2 

1=2  VEGETATION  LAYER  TWO  J=3 

1=3  VEGETATION  LAYER  THREE  J=4 

GROUND  LAYER  J=5 


Figure  1.  Diagram  showing  sink  and  source  variable  indices  used  in 

in  the  model  energy  flow  formulations 
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Figure  8.  Concluded. 
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Figure  10.  Leaf  area  Index  distribution  for  the  Douglas-fir  canopy 
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11.  Leaf  area  index  distribution  for  the  oak-hickory  canopy 
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Figure  12.  Measured  average  layer  needle  temperatures 
for  Douglas-fir  plotted  with  air  and  ground  temperatures 
for  a  48-hour  period  from  August  4  to  August  5,  1979; 
layer  1  and  layer  2  needle  temperatures  measured  at  26  m 
and  20  m  respectively 
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Figure  13.  Measured  average  canopy  temperature  for  oak- 
hickory  plotted  with  air  and  ground  temperatures  for  a 
48-hour  period  from  August  18  to  August  19-  1979;  canopy 
temperatures  measured  intermittently  from  1100  hours  on 
August  18  to  1200  hours  on  August  19  with  a  hand-held 
thermal  radiometer 


t 


(vW/n)  aiis 


Figure  14.  Plots  of  thermal  model  environmental  input 
parameters  for  the  Cedar  River  site  from  0000  hr  4  August 
1979  to  2400  hr  5  August  1979:  a)  Global  shortwave  radi¬ 
ation  (SWR),  b)  Wind  speed  (WS),  c)  Relative  humidity  as 
estimated  from  wet  and  dry  bulb  temperatures,  and  d)  air 
temperature  (  Q  )  and  ground  temperature  (  O  ) 


Figure  15.  Plots  of  thermal  model  environmental  Input 
parameters  for  Walker  Branch  Site:  a)  Global  shortwave 
radiation  (SWR),  b)  Wind  speed  (WS),  c)  Relative  humidity 
as  estimated  from  wet  and  dry  bulb  temperatures,  and 
d)  air  temperature  {  Q  )  and  ground  temperature  (  O  ) 
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Figure  18.  Layer  2  predicted  temperature  plotted  with 
average  temperatures  measured  at  the  20-m  level  In  the 
Douglas -fir  canopy 
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Figure  19.  Simulation  results  for  the  three-layer 
oak-hickory  canopy  plotted  with  air  temperatures 
for  the  48-hour  period  from  18-19  August  1979 
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Figure  20.  Layer  2  predicted  temperature  plotted  with 
measured  average  temperature  of  oak-hickory  canopy; 
measurements  made  with  a  thermal  infrared  radiometer 


APPENDIX  A:  PROGRAM  LISTINGS 


TMODEL 


PROGRAM  TMODEL 1 ( INPUT , OUTPUT ,  TAPE3*INPUT , TAPE4*0UTPUT  , 
1TAPE1 ,  TAPE2,TAPE3) 


t*********THODEL--VERSION  2. 14. BO************** 
COMMON/PARA1/SIG,STA(3) , S(3,3) ,ST6(3) ,X(3),A(3) ,U(4,9) » 
1ALP<3),FMU,TA,TAC,TG,EPS<3),EPSTG,RM,RL 
C0MM0N/PARA2/BTA, BTG,BX<3) , DBX(3) f QX<3) ,DQX<3) »RX(3),DRX<3) 
COMMOM/EST IM/FX < 3 ) #  DFX ( 3 , 3 ) 

C0NM0N/SENS0R/ERT<9) , ERTH(3) ,EX<9) ,EXH<3> 

DIMENSION  DX<3) ,SABS(3) 


SIG*3.6686E-8 


TOL*. OOOOI 


READ  CANOPY  GEOMETRY  MATRICES 
FROM  TAPE  2 


READ(2,203)  TITLE1 ,TITLE2, TITLES 
URITE  <  3 , 203 ) TITLE 1 tTITLE2f TITLE3 
C  READ  U  MATRIX:  SENSOR  VIEW  ANGLE  UEIGHTS 

C 

DO  2001  N«1,4 
READ(2,202)  <IKN, 

2001  CONTINUE 
C 

C  —  READ  S  MATRIX:  LU  FLUX  TRANSFERS 
C 

DO  2000  1*1 f 3 

READ<2f202)STA(I), <S<It J> ,J*1 ,3) fSTG(I) 

2000  CONTINUE 
C 

C  —  READ  ABS:  Stf  FLUX 
C 

202  F0RNAT<9F7.4> 

RE AD (2, 202)  (SABS(N) fN*1 ,3) 

203  F0RMAT(1Xf3A10) 

C  —  DEFINE  NOMINAL  VALUES  FOR  MODEL  PARAMAMETERS 
C 

C1000  FORMAT ( 14 f4F5.0,2XtF5.0) 

C1000  F0RHAT(I4,2F7.2f 2F0.2, 2XfF10.2> 

C 

ALP(1)*1.0 
ALP(2)*1 .0 
ALP(3)*1 .0 
EPS<1)*1.0 
EPS(2)*1 .0 
EPS(3)*I.O 
EPST6-1.0 


/ 


CJ  o 


—  READ-IN  THE  NUMBER  OF  SIMULATION  PERIODS  — 

PRINT  201 

201  F0RMAT<//,*  ENTER  THE  NUMBER  OF  SIMULATION  PERIODS  DESIRED*/) 
READ*, NSIM 
C  —  READ  IN  CANOPY  RESISTANCE  — 

PRINT  210 

210  FORNAT(/,*EKTER  THE  CANOPY  STOMATAL  RESISTANCE  FOR  THIS  RUN*/) 
READ*,RL 
URITE(6,600) 

000  FORMAT ( 1H  ,4HTINE,4X,3HSWR,15X,1HA,31X,1HB,31X,1HH,23X,2HLE) 
WRITE (A,A02) 

002  FORMAT (1H  ,17X,1H1 ,7X,1H2,7X,1H3,7X,1HA,7X, 1H1 ,7X,1H2,7Xf 1H3,7X, 
.1H6,7X, 1H1 ,7X,1H2,7X,1H3,7X, 1H1 ,7X, 1H2,7Xf 1H3) 

2  CONTINUE 

DO  100  NTIHE*1 ,NSIN 

READ ( 1 ,*)ITIME,TA,T6f FMU,RH,GLB 

C .  CONVERT  FMU  IN  H/SEC  TO  CM/SEC . 

FMU  «  FHU*100. 

A ( 1 )=GLB*SABS( 1 )  $  A<2)«6LB*SABS<2)  t  A(3)  =  6LB*SABS(3) 

TAC-TA 
DO  00  1*1 ,3 
90  X(I)*TA 
50  CALL  FEVAL 
DO  20  1*1,3 
20  FX(I)*-FX(J) 

CALL  SOLVE(FX,3,DFX,3,DX) 

DO  30  1*1,3 
30  X(I)=X(I)+DX(I) 

DO  40  1*1,3 
DEV*DX(I) 

IF( ABS(DEV)  .6T.  TOL)  GO  TO  50 
40  CONTINUE 
CALL  MATTS 

WRITE ( 3 , 80 >  IT I ME , T A , TG , (ERTH< J),J=1 ,3) 

BO  FORMAT ( 1H  r 1 1 0 , 2F7 . 2 , 3F 6 - 1 ) 

05  FORMAT (I10,4F10.5> 

BTA*BTA«SIG  t  BX ( 1 )=BX( 1 ) *SIG  ♦BX(2)*BX(2)*SIG  4  BX(3)*BX(3)*SI6 
BTG*BT0*SI6 
C 

MRITE(6,604)ITIME,GLB, ( A( I ) , 1  =  1 ,3) ,BTA, (BX( I ) ,  1*1 , 3) , BTG , (QX ( I ) , 
.1*1 ,3) , (RX< 1) , 1*1 ,3) 

404  FORMATdH  ,14, 15F8.2) 

100  CONTINUE 
STOP 
END 

SUDROUTINE  FEVAL 
C 

C0MM0N/PARA1/SIG,STA(3) ,S(3,3),ST6(3),X(3),A(3),U(4,0), 

1 ALP<3) ,FNU,TA,TAC,T0,EPS(3) ,EPSTG,RH,RL 
C0MH0N/PARA2/BTA,BTG,BX<3) ,DBX(3) ,QX(3) ,DQX(3> ,RX(3) , DRX ( 3 ) 


c 


C0HH0N/ESTI!1/FX<3),BFX<3,3> 


CALL  FBTA(TA,3TA) 

CALL  BFUNC(cP3TG,TG,2TG,D3TG> 

DO  10  1*1,3 

CALL  3FUNC(EPS(  I ),X(I),3X(I),D8X(D) 

CALL  QFUNC(X( I)  ,TAC,F!1U,QX( I )  ,DQX ( I ) ) 

CALL  RFUNC ( X( I ) , TAC ,f HU, SL , RH , RX ( I ) , DRX ( 1 ) ) 

10  CONTINUE 
DO  20  IL*1 ,3 

20  FX<IL)=0.5*ALP(IL)*SIG*<BTA*STA<IL)*3X<1 )*S<  IL,  1  )+3X( 2) * 
1S(IL,2)+BX(3)*S(IL,3)+BTG*STG(IL))+A(IL)-SIG*3XUL)+3X(IL) 
2+RX1IL) 

DO  30  1*1,3 
DO  30  J*1 ,3 
30  DFX<I,J>*0. 

DO  40  I L * 1 ,3 
DO  40  J*1 ,3 
IF(J.NE.IL)  GO  TO  35 

DFX(IL,J)*0.3*AL?<IL>*SIG*DBX<J)*S<IL,J)-5IG*BBX(IL)+DQX<IL)* 
1 DRX ( I L ) 

GO  TO  40 

33  DFX(IL,J)*0.5«ALP<IL)*SIG*D8X( J)*S(IL,J) 

40  CONTINUE 
RETURN 
END 

SUBROUTINE  3FUNC<£?SI, XI, 3X1,33X1) 

3XI*EPSI» (XI *273.0) **4. 

DBXI=4.*EPSI*(XI*273.0)**3. 

RETURN 

END 

SUBROUTINE  F8TA(TA,3TA) 

E?STA*1.-0.241*EXP(-7.77E-4*TA*TA) 

CALL  BFUNCtEPSTA, TA,8TA,DBTA) 

RETURN 

END 

SUBROUTINE  QFUHCUI  ,TAC,FNU,OXI  ,30X1 ) 

IF (FHU.ST.30. )  GO  TO  10 
HC=0.4?775*<20.4*Q,2*FNU**0.97) 

GO  TO  20 

10  HC*0.6??75*(0.?5*FNU**0.97) 

20  QXI=(XI-TAC)*(-HC> 
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08XI*-HC 
RETURN 
END 

SUBROUT I NE  RFUNC <  X I , TAC , FNU ,  RL , RH , RX I ,  DRX I ) 

RNUN«FEX(XI>*1 .0E-6-RH*FEX<TAC)*1 .OE-6 
RDEN*RL+( 1 ./60. )*<  0,04+1 .27*FHU**<-0.5> ) 

RX1 1 »-697.?5*< -0.566*X1*5?7.3) 

RXI2*RNUH/RDEN 
RX I *RX 1 1  *  R  X 1 2 

DRXI*6?7.75*0. 566*RXI2+RXI1 *<0.05671 5E-6*FEXai)>/RBEN 
RETURN 
END 

SUBROUTINE  INVERSE( AtN,D> 

—  INVERT  A  3*3  REAL  MATRIX  A  UHQSE  DETERMINANT  IS  D 
THE  RESULT  UILL  BE  STORED  IN  A 

DIMENSION  A<3,3),3<3,3> 

C 

D*A<1,1)*A<2,2)*A<3,3)+A<1,2>*A<2,3>*A<3fl>+AU,3)*A<2,1>* 
IA<3,2)-A<3f1)*A<2,2)*A<l,3)-A<1,1>*A(3,2)*A(2f3>-A<2,1)* 
2A( t ,2)*A(3,3) 

B<1,1)*<A(2,2)*A<3f3)-A<2f3)*A<3,2)>/B 
B< 1 ,2)*-< A< 2, 1 )*A<3,3)-A(2, 3)*A< 3, 1 ) >/D 
B<1,3)*<A<2,1)*A<3,2)-A<2,2>»A<3,1>>/D 
B(2,1)*-<A(1,2)*A<3,3)-A<I,3)*A(3,2))/D 
B(2,2)»(A<1f1)*A(3t3)-A(1,3)*A<3,1))/D 
B(2,3)*-(A(1,1 )*A<3r2) -A( 1 f 2)*A(3, 1 > )/D 
B<3,1)=(A(1,2)*A<2,3)-A<1f3)*A<2,2>)/D 
B<3,2)—<A<1,1)*A<2,3)-A<1f3)*A<2,l))/D 
B<3,3)*<A<1,1)*A(2,2)-A(l,2)*A(2f1))/D 
DO  10  1*1,3 
DO  10  J*1 ,3 
10  A<I,J)*B<I,J) 

RETURN 

END 

FUNCTION  FEX(XI) 

C 

XX*5. 2342*EXP<  0.05671 5*X I ) 

FEX*XX 

RETURN 

END 

SUBROUTINE  SOLVE< T,N, A,NfX) 
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C  DIMENSION  Y(N),A<N,N),X<M),ATA<N,N),ATY(N) 
DIMENSION  Y(3),A(3,3),X<3),ATA(3,3),ATY(3) 
C 

DO  tO  1*1  ,N 
DO  tO  J*1 ,N 
ATA<I, J)*0. 

DO  tO  K*t,N 

10  ATA<I, J)*ATA<I,J)+A<K,I)*A<K,J> 

CALL  INVERSE* ATA, N,D> 

DO  20  1*1,  M 
ATY  < I )*0. 

DO  20  J*t,N 

20  ATY(I)*ATY(I)+A(J,I>*Y<J) 

DO  30  1*1, M 
X(I)*0. 

DO  30  J=1,M 

30  X<  I  )*XU  )+ATA<  I ,  J>*ATY <  J) 

RETURN 

END 


SUBROUTINE  UATTS 

COMMON/PARAI/SIG,STA( 3) ,S(3,3> , 3TG( 3) , X ( 3) , A ( 
1  ALP  <  3 ) ,FNU,TA,TAC,TG,EPS<3),EPSTG,RH,RL 
C0MH0N/SENS0R/ERK9)  ,ERTH(3)  ,EX( 9) ,EXH<  3) 

C 

DO  1  N*1 ,3 

EXH<N)»EPS(M)*SIG*CX(N) +273. 1**4 
ERTH(N)»<  (EXH(N)/SIG)**0.25)  -  273. 

1  CONTINUE 
C 

DO  2  J»1 ,9 

EX(J)*SIG*EPS(1)»Wn,J)*<X(1  )+273.)**4  + 
1SIG*EPS<2)*U(2,J)*<X(2)+273.)*»4+ 
2SIG*EPS<3)*U<3,J>*<X(3H-273.t**4+ 
3SIG*EPSTG*U(4, J)*<TB*273. )**4 
ERT < J ) « <  <EX(J)/SIG>**0. 251-273. 

2  CONTINUE 
RETURN 
END 
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SCALC 


PROGRAM  SCALC< INPUT, OUTPUT, TAPE5=INPUT,T APE6=OUTPUT, TAPE  1 ,TAPE2, 
.TAPE3) 

r* 

l 

C*»*»***#****»»SCALC--VERSI0N  2.14.80  ***************** 

C 

C 

C . . .THIS  PROGRAM  IS  EXTRACTED  FROM  THE  TCSM  MODEL  BY  KIHES  TO 
C... CALCULATE  THE  SENSOR  V1EV  WEIGHTING  MATRIX  U  AND  THE  LONGWAVE 

C _ RADIATION  EXCHANGE  MATRIX  S  FOR  A  GIVEN  CANOPY  GEOMETRY. 

C.. .REQUIRED  INPUTS  ARE  LEAF  INCLINATION  ANGLE  AND  CORRESPONDING 
C.. .FREQUENCY,  CANOPY  DENSITY  < SLAI )  AND  LEAF  AREA  INDEX  (FLAI ) 

C...FOR  EACH  OF  THREE  CANOPY  LAYERS. 

C 

C 

COMMON  /A/  UV,RH,RL,D(3) 

COMMON  HI  PGAP<3,9),  PHIT«3,9)  ,PGAP2<3,9) ,  PHIT2<3,9> , STEF 
COMMON/C /COST A (9,9,18),EHISSV(4),ABS0RB(3),  ESKY,  EGRD,  SECTAR(9) 
COMMON/ D/CONT <Z,5,9),C(3,5,9),  SUMT(3,9),  KELV,  GT,  NUSIM,  ITIME 
COMMON/E/AT, THETA < 9) , PHI ( t 8 ) ,XLF ( 9 ) ,  YLF (9) ,  ZLF(9),  XS(9,18) 
C0MM0N/F/YS(9,18),2S(9),  CEDTR,  B,  FREQD(9,3)  ,  UA( 15) ,EPS 
COMMON/G/NSIG,N,  ITMAX 
COMMON  /H/  INDEX1 ,TITLE<85 
COMMON  III  X < 3 ) 

COMNON  /J/  THERM, THMLEX ,CONVEC, TRANS 
COMMON/S/  ABS0L<3) 

C 

C...READ  AND  ASSIGN  THE  INPUT  DATA 
C 

NUSIM=-1 
INDEX1 »0 

76  IF ( INDEX  1 . EO . NUSIM ) STOP 
CALL  INPUTDA 

I F « INDEX  1 . GT . 1 )  GO  TO  95 
C 

C... CALCULATE  THE  CANOPY  GEOMETRY  COEFFICIENTS 
C 

CALL  CANGEOM 
C 

CALL  WHAT (TITLE) 

C 

C... CALCULATE  THE  SIN  THETA  FACTORS  FOR  ALL  SOURCE  ANGLE-LEAF  ANGLE 
C  PURMUTATIONS. 

C 

CALL  DEVANG 
C 
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C... CALCULATE  THE  NORMALIZING  FACTOR  FOR  THE  RELATIVE  SIZES  OF  SOURCE  SECTORS 
C 

CALL  SECTOR 
95  CONTINUE 
C 

C.. .CALCULATE  THE  THERHAL  RADIATION  COEFFICIENTS 
C 

CALL  SETUP 
C 

C... CALCULATE  THE  AVERAGE  LEAF  TEMPERATURE  WITHIN  EACH  LAYER. 

C 

C... DISPLAY  THE  OUTPUT 
C 

CALL  OUTDAT 
GO  TO  76 
END 

SUBROUTINE  INPUTDA 
C 
C 

C... SUBROUTINE  INPUTDA  READS  AND  ASSIGNS  THE  INPUT  DATA 

C 

C 

COHHON/GEO/  PHIT 1 ( 3 , 9 ) , FLAI ( 3, 1 ) ,SLAI ( 3, 1 ) , AXLFA ( 1 9 . 3 ) , AYLFA < 1 V , 3 ) 
COMMON  /A/  UV , RH , RL , B ( 3  ) 

COMMON  III  PGAP<3,9),  PHIT<3,9) ,PGAP2<3,9) ,  PHI T2 ( 3 , 9 ) , STEF 
C0MM0N/C/C0STA<9, 9, 18), £MISSV<4>, ABSORBS),  ESKY,  EGRL ,  SECTARY} 
COMMON/D/CONT (3,5,9),C(3»5,9),  SUMT<  3,9) ,  KELL*,  GT ,  NUSIM,  I T I  ME 
COMMON/E/AT, THETA19 > ,PHI < 1 8> ,XLF (9) ,  YLF ( 9  > ,  ZLF(9>,  XS(9,18) 
C0HM0N/F/YS(9,18),ZS<9),  CEBTR ,  B,  FREQB(9,3>  ,  UA ( J 5 ) , EPS 
CGMMON/G/NSIG , N ,  ITMAX 
COMMON  tW/  INDEXt , TITLE (E) 

COMMON  III  X(3) 

COMMON  m  THERM, THMLEX , CONVEC , TRANS 
COMMON  /HI  ST0R(3) 

COMMON/S/  ABSQL(3) 

C 

C. . .TEST  FOR  THE  SIMULATION  NUMBER  AND  SKIP  TO  THE  APPROPRIATE  INPUT  DATA 
C 

IF  (INDEX!. EO.O)  GO  TO  99 
IF(INDEX1.EQ. NUSIM)  STOP 
99  CONTINUE 
C 

C.. .ASSIGN  THE  STEFFAN  BOLTZMANN  CONSTANT  UATTS/M**2*K**4 
C 

STEF*5.6686E-8 
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C... ASSIGN  THE  CONVERSION  FACTOR  FOR  KELVIN-DEGREES 
C 

B*  273.0 
C 

C...READ  THE  AVERAGE  THERMAL  EM1SSIVITY  COEFFICIENTS  FOR  THE  3  VEGETAION  LAYER 
C  (1,2,3)  AND  THE  GROUNDH). 

C 

READ( 1 , 199) (TITLE(N) ,N*1 ,8) 

IF  <  EOF ( 1 ) .NE.O. ) S  TOP 
199  FORMAT (6A10) 

C 

C 

C 

C 

C. . .READ  THE  CANOPY  GEOMETRY  FREQUENCY  DISTRIBUTIONS  OF  THE  ELEMENTS 
C  IN  LAYERS  1,2,3.  AXLFA  REPRESENTS  THE  INCLINATION  ANGLES  0-90 
C  (5  DEGREE  INTERVALS)  AND  AYLFA  REPRESENTS  THE  CORRESPONDING 
C  FREQUENCY.  SLAI  AND  FLAI  ARE  EACH  LAYERS  S  PARAMETER  AND  LAI 
C  RESPECTIVELY. 

C 

C 

DO  190  1*1,3 

READd ,*)(AXLFA(M,I),AYLFA(M,I),M*1 ,19) 

READ ( 1 ,*)SLAI (1,1) , FLAK  1 ,1 ) 

190  CONTINUE 

101  FORMAT  (8F10.5) 

RETURN 

END 

SUBROUTINE  OUTDAT 
C 

c 

C... SUBROUTINE  OUTPUT  FORMATS  THE  DATA  TO  BE  DISPLAYED. 

C 

C 

COMMON/SENS/  ELA YT  <  9 ) ,ELA YH ( 3 ) , ERTT ( 9 ) , ERTH ( 3 ) 

COMMON/GEO/  PHIT1 (3 , 9) , FLAI (3, 1 ) ,SLAI ( 3, 1 ) ,AXLFA( 1 9,3 ) , A YLFA( 1 9, 3 ) 

COMMON  /A/  UV,RH,RL,D(3) 

COMMON  III  PGAP<3,9),  PHIT (3, 9)  ,PGAP2(3,9) ,  PHIT2(3,9) ,STEF 
COMMON/C/COST A < 9, 9, 18),EHISSV(4),ABS0RB(3),  ESKY,  EGRD,  SECTAR ( 9 ) 
CQHHON/B/CONT (3,5,9),C(3,5,9>,  SUMT(3,9),  KELV,  GT,  NUS1M,  ITIME 
COMMON/E/AT, THETA ( 9) ,PHI < 1 8 ) ,XLF (9) ,  YLF ( 9) ,  2LF(9),  X S < 9 , 1 8 ) 
C0MM0N/F/YS(9,18),ZS(9),  CEDTR,  B,  FREQH (9,3)  ,  UA(15),EPS 
COMMON/G/NSIG,N,  ITMAX 
COMMON  /H/  INDEX1 ,TITLE<6) 
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COMMON  /I/  X ( 3 ) 

COMMON  / J/  THERM, THMLEX, CON VEC, TRANS 
COMMON  /K/  TT1(3),TT2<3),TT3<3),TT4(3) 

COMMON  /N/  STOR(3) 

COMMON  /L/  TEMP ( 3 ) 

COMMON/S/  ABSOL ( 3 ) 

DIMENSION  S ( 3 , 5 ) 

C 

c 

C.. .WRITE  THE  CALCULATED  GEOMETRY  FOR  EACH  LAYER 
C 

DO  319  1=1,3 
WRITE(2,320)  I 

320  FORMAT  (///,*  THE  COMPONENT  ANGLE  COMPUTATIONS  FOR  LAYER  * . 1 1 , / ) 
URITE(2,321 )  FLAI (1,1), SLAI (1,1) 

321  FORMAT  (  *  LAI  =  *,F4.2,4X,*  S=  *,F4.2,/> 

WRITE (2, 322)  (AXLFA(M,I),AYLFA(M,I),M=1,19) 

322  FORMAT (  *  XLFA,YLFA  *,/,  (2X,UF8.3) ) 

UF:ITE(2,323)  (PGAP< I ,M>  ,M=1 ,9) 

323  FORMAT <//,*  PGAP  FOR  1-9  INCLINATION  INTERVALS*, 9F8. 3) 

319  CONTINUE 

C 

C... WRITE  THE  CALCULATED  THERMAL  CONTRIBUTIONS  COEFFICIENTS 
C 

WRITE<2,302) 

302  FORMAT (IX,/,*  THE  PROPORTION  OF  RADIANCE  AREA  CONTRIBUTED  BY 
«A  SECTOR  OF  THE  9  BANDSU-9)  DIVIDED  BY  18  (SECTORS)  ARE**,/) 

WRITE (2,303 ) (SECT ARC  I ) , 1  =  1 ,9) 

303  FOF:MAT(10X,?F10.5,//) 

WRITE(2,40) 

40  FORMAT  (IX,//,*  THE  BAND-PGAP-PHIT-COEr FICIENTS  FOR  THE  THERMAL  F.'A 
♦BIATION  TRANSFERS  ARE  =*,/) 

DO  39  1=1,3 
URITEC2.41)  I 

41  FORMAT  (IX.*  THE  9  BAND  COEFFICIENTS  TO  LAYER  *,I1,*  ARE*) 

DO  39  J=1 ,5 

WRITE  (2,42)  J,  (CONT ( I , J,M) ,M=1 ,9) 

42  F0RMAT(8X,*  FROM  L AYER* , 1 1 ,2X , 9F6. 4 ) 

39  CONTINUE 

WRITE(2,50) 

50  FORMAT  (IX,////,*  THE  FINAL  THERMAL  RADIATION  COEFFICIENTS  ARE  AS  FOLLOWS 

+*,/) 

DO  49  MXX=1 ,3 
DO  49  NXX  =  1 ,5 
S(MXX,NXX)  =  0.0 
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49  CONTINUE 
DO  51  1=1,3 
URITE<2,52)  I 

52  FORMAT  (IX,*  THE  THERMAL  RADIATION  CONTRIBUTION  TO  LAYER  *,I1,*  FO 
+R  EACH  OF  THE  9  LEAF  INCLINATIONS  ARE*) 

DO  51  J=1 ,5 

WRITE  (2,53)J,(C(I,J,M),M=1,9) 

DO  51  K=1 ,9 

S(I,J)  =  S(I,J)+C(I,J,K)*FREQD(K,I) 

53  FORMAT  (8X,»  FROM  LAYER*, II , 2X, 9E1 0.3 ) 

51  CONTINUE 

DC  55  IXX  =  1,3 

URITEl 3,505 )(S(IXX,JXX),JXX=1 ,5) 

505  FORMAT (5F7.4 ) 

URITE(6,503)(S(IXX,JXX),JXX=1 ,5) 

503  FORMAT ( 1 H  ,5F10.4) 

55  CONTINUE 
RETURN 
END 

SUBROUTINE  SETUP 
C 
C 

C... SUBROUTINE  SETUP  PRE-CALCULATES  AND  PRE-ARRANGES  MANY  OF  THE  THERMAL 
C  COEFFICIENTS  NEEDED  FOR  THE  FINAL  ENERGY  BUDGETS  UHICH  ARE  PLACED  INTO  THE 
C  2SYSTM  ROUTINE. 

C 

C 

COMMON  /A/  UV,RH,RL,D(3) 

COMMON  III  PGAP (3,9),  PHIT ( 3, 9) , PGAP2<3, 9 ) ,  PHIT2 ( 3,9) , STEF 
COHMQN/C/COSTA (9, 9, 1 8 ) , EMISSV (4 ) , ABSORB (3 ) ,  ESKY,  EGPD,  SECTARY) 
C0MM0N/D/C0NT(3,5,9),C(3,5,9),  SUKT(3,9),  KELV,  GT,  NUSIM,  I T I ME 
COMMON /E/AT, THETA (9 > ,PHI(1B),XLF(9),  YLF ( 9) ,  ZLF<9),  XS(9,1B) 
C0HM0N/F/YS<9,18),ZS<9),  CEDTR,  B,  FREGB(9,3)  ,  UA(15),EPS 
COMMON/G/NSIG , N ,  ITMAX 
COMMON  /H /  INDEX1 , TITLE < 8 ) 

C 

C 

C...FOR  EACH  LAYER  CALCULATE  THE  BAND-PGAP-PHIT  COEFFICIENTS  NEEDED  FOR  EACH 
C  LAYERS  THERMAL  RADIATION  CONTRIBUTION  TO  A  SPECIFIC  LAYER. 

C 

C 

DO  20  1=1,9 
C 
C 

C.. .CONTRIBUTION  COEFFICIENTS  TO  LAYER  1 


All 


( 


c 

c 

c.. 

c 

c 

c.. 

c 

c 

c.. 

c 

c 

c.. 

c 

c 

c.. 

c 

c 

c 

c.. 

c 

c 

c.. 

c 

c 

c.. 

c 

c 

c.. 

c 

c 

c.. 

c 

c 

c.. 

c 


...FROM  SKY 

CONT(1,1,I)=  PGAP2<1,I> 

...FROM  LAYER  1 

CONT (1,2,1)=  2.*PHIT2< 1 , 1 ) 

...FROM  LAYER  2 

CONT (1,3,1)=  PGAP2 ( 1 ,1)-PGAP2(1 , 1 )*PGAP( 2, 1 ) 

...FROM  LAYER  3 

CONT< 1,4,1)=  P6AP2( 1 , 1 )*PGAP(2, 1 )-PGAP2< 1 , 1 ) *PGAP(2, 1 )*PGAP<3, 1 ) 
...FROM  GROUND 

CONT < 1 ,5,1 >=  PGAP2U , I >*PGAP< 2 , 1 )*PGAP< 3, 1 > 

•CONTRIBUTION  COEFFICIENTS  TO  LAYER  2 

...FROM  SKY 

CONT (2,1,1)=  PGAP<  1 ,1  )*P6AF'2(2, 1 ) 

...FROM  LAYER  1 

CONT (2,2,1)=  PGAP2(2 , I )-PGAP2( 2, I )*PGAP( 1,1) 

...FROM  LAYER  2 

CONT (2,3,1 )=  2. *PHIT2<2, 1 ) 

...FROM  LAYER  3 

CONT (2,4,1)*  PGAP2(2,I)-PGAP2(2,I)*P6AP<3,I> 

...FROM  GROUND 

CONT (2,5,1)=  PGAP2 ( 2, 1 )*PGAP< 3, 1 ) 


A12 


T 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


. . .CONTRIBUTION  COEFFICIENTS  TO  LAYER  3 


...FROM  SKY 

CONT (3,1,1)* 

...FROM  LAYER  1 

CONT  (3,2, 1 )« 

...FROM  LAYER  2 

CONT (3,3,1)= 

...FROM  LAYER  3 

CONT (3,4,1)* 

...FROM  GROUND 

CONT (3,5,1  )  = 
20  CONTINUE 


PGAP< 1 , 1 )*PGAP(2, I )*PGAP2(3, I ) 


PGAP2!3,I)*PGAP(2,I)-P6AP2<3,I)*PGAP(2,I)»PGAP(1 


PGAP2(3,I>-PGAP2<3,I>*PGAP(2,I) 


2.*PHIT2(3, I ) 


PGAP2 (3,1 ) 


I) 


C 

C 

C...NOU  FORK  THE  EQUATION  COEFFICIENTS  FOR  THE  CONTRIBUTED  THERMAL  RADIANT 
C  ENEP.'GY  TO  EACH  LAYER  AND  FOR  EACH  LEAF  INCLINATION  ANGLE  UITHIN  A  LAYER. 
C 
C 

CALL  SET03<C, 3,5,9) 

C 

C...THERHAL  RADIATION  CONTRIBUTION  TO  LAYER  N 

r 

w 

DO  30  N=1 ,3 
C 

C...FOR  EACH  LEAF  INCLINATION  ANGLE  INTERVAL 
C 

DO  30  I«  1,9 
C 

C...SUK  EACH  SECTORS  RADIATION  CONTRIBUTION  (9  BANDS  CONTAINING  18  SECTORS) 

C 

DO  30  J=1 ,9 
DO  30  K=1 , 1 6 
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c 

C... ABSORBED  THERHAL  RADIATION  CONTRIBUTED  BY  SKY 
C 

CYN, 1 ,1 )*  C  ( N ,  1 , 1 )  +  SECTARY  J)*COHT <H, 1 ,J) 

+*COSTA( I, J,K) 

C 

C.. .ABSORBED  THERHAL  RADIATION  CONTRIBUTED  BY  LAYER  I 
C 

C(N,2,I>=  C(N,2, I )  ♦  SECTAR(J)*C0NT(N,2,J) 

+*CQSTA(I, J,K) 

C 

C... ABSORBED  THERHAL  RADIATION  CONTRIBUTED  BY  LAYER  2 
C 

C(N,3,I)=  C(N,3, 1 )  +  SECTAR(J)*C0NT(N,3,J) 

+*COSTA ( I ,  J ,  K ) 

C 

C.. .ABSORBED  THERHAL  RADIATION  CONTRIBUTED  BY  LAYER  3 
C 

C(K,A, 1)=  C(N,4, 1)  ♦  SECTARf J)*CONT(N,4,J) 

♦*COSTA(I, J,K) 

C 

C... ABSORBED  THERHAL  RADIATION  CONTRIBUTED  BY  THE  GROUND 
C 

C(N,5,I)=  C(N,5,n  ♦  SECTARY J)*CONT (N,5, J) 

♦*COSTA ( I , J,K ) 

30  CONTINUE 
RETURN 
END 

SUBROUTINE  DEVAN6 

r 

u 

C.. .SUBROUTINE  DVANG  CALCULATES  THE  COS(ANGLE)  DEVIATION 

C  ANGLE  OF  ALL  LEAF  INCLINATIONS  SOURCE  0RIENTAI0N5  PERHUTATIONS.  THE  THEOK 
C  IS  BASED  ON  THE  EXISTENCE  OF  PLANE  ELEHENTS  AS  USED  IN  THE  SRVC  HODEL. 

C 

C 

COHHON  /A/  UV,RH,RL,D(3> 

COHHON  /B/  PGAP  Y  3  f  9 ) ,  PHIT ( 3 , 9 ) , PGAP2 < 3 , 9 ) ,  PKIT2 ( 3 , 9 )  r  STEF 
COHHQN/C/CQSTA (9,9,18) ,EHISSV(4) f ABSORB ( 3 ) ,  ESKY,  EGRD,  SECTAR Y 9 ) 
COHHON/D/CONT (3, 5,9) ,C(3,5,9) ,  SUHT(3,9),  KELV,  GT,  NUSIH,  ITIHE 
COHHON /E/AT, T HE TA(9),PHI(18>,XLF(9),  YLF ( 9) ,  2LF(9),  XS  <  9 , 1 8 ) 
COHHON/F/YS( 9, 1 8 > ,ZS< 9 ) ,  CEDTR,  B,  FREQD (9,3)  ,  «A Y 1 5 ) , EPS 
COHHON/G/NSIG,N,  ITHAX 
INTEGER  SB,SS 
CEDTR*  0.01 7453293 


c 

C... CALCULATE  INCLINATIN  ANGLES  IN  RADIANS 
C 

THETA ( 1 )  =  5.  *  CEDTR 
DC  fO  1=1,8 

THETA(I  +  1 )=  T HE T A ( I )  ♦  10.0  *  CEDTR 
10  CONTINUE 
C 

C.. .CALCULATE  AZIMUTH  ANGLES  IN  RADIANS 
C 

PHI ( 1 )=  10.*CEDTR 
DO  20  1*1 ,17 

PHI (1  +  1  )=  20.*CEDTR+PHI(I) 

20  CONTINUE 
C 

C... CALCULATE  ALL  THE  DIRECTION  COSINES  OF  SOURCE  SECTORS 
C 

DO  AO  1*1,9 

ZS< I )=SIN( THETA(  I ) ) 

BO  AO  J=1 , 1 8 

XS<I,J)=COS(THETA(I ) )*COS(PHI ( J) > 

YS(!,J)=  COS(THETA(I))*SIN(PHI(J) ) 

AO  CONTINUE 
C 

C.. .CALCULATE  THE  DIRECTION  COSINES  FOR  THE  NORMAL  VECTOR  OF  ALL  PLANAR  LEAF 
C  INCLINATION  ANGLES  ASSUMING  THAT  THE  AZIMUTH  ANGLE  IS  EQUAL  TO  ZERO  DEGREES . 
C 

DO  30  1=  1,9 

XLF ( I  )  =  -SIN ( THETA ( I ) ) 

YLF( I )=  0.0 

ZLF  <  I  )=  COS  <  THETA  < I ) ) 

30  CONTINUE 
C 

C... CALCULATE  THE  ABSOLUTE  VALUE  OF  THE  DOT  PRODUCTS  OF  ALL  SOURCE-LEAF 
C  ANGLE  PERMUTATIONS.  THIS  VALUE  IS  EQUAL  TO  THE  COSINE  FACTOR  DESIRED. 

C 

DO  50  11*1,9 
DO  50  SIM, 9 
DO  50  SS=  1,18 

DOT*  <XLF<LI)*XS<SB,SS)+YLF<LI)*YS(SB,SS)+ZLF(LI)*ZS(SB>> 

COSTA(LI ,SB,SS)=  ABS  (DOT) 

50  CONTINUE 
RETURN 
END 

SUBROUTINE  CANGEOM 
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c 

c 

C... SUBROUTINE  CANGEQM  CALCULATES  THE  CANOPY  6E0HETRY  COEFFICIENTS. 

C...THE  SUBROUTINE  CANGEON  CALLS  SUBROUTINE  SRVCMOD  UHICH  IS  A  MODIFIED 
C  PORTION  OF  THE  SRVC  HODEL  THAT  CALCULATES  THE  CANOPY  GEOMETRY 
C  PARAMETERS. 

C 

C 

COMMON/GEO/  PHIT1 <3,9 ) ,FLAI <3, 1 ) ,SLAI ( 3, 1 ) f AXLFA (1 9,3) , AYLFA< 1 9 ,3 ) 
COMMON  /A/  UV , RH , RL , D ( 3 ) 

COMMON  /B /  PGAP(3, 9) ,  PH1T ( 3 , 9 > ,PGAP2<3,9>,  PHIT2(3,9) f STEF 
C0NN0N/C/C0STAY9, 9, 18), EMISSV<4), ABSORBS),  ESKY,  EGRD ,  SECTAR(9> 

COMMON /Ii/CONT  <3,5,9),C<3,5,9),  SUMT(3,9),  KELM,  GT,  NUSIM,  I T I  ME 
C0HM0N/E/AT,THETA<9),PHI<18>,XLF(9>,  YLF (9) ,  2LF  <  9 ) ,  XS<9,18) 
COMhOH/F/YS < 9,18)tZS<9),  CEDTR,  B,  FREQD<?,3)  ,  UA(15),EPS 
COHMQN/G/NS IG , N ,  ITMAX 
CALL  SRVCMOO 
DO  10  1=1,3 
DO  10  M=1 , 9 
C 

C.. .TRANSFER  IDENTICAL  ARRAYS  PHIT  AND  PHIT1 .  PHIT  CONTAINS  THE 
C  PROBABILITY  OF  HIT  COEFFICIENTS  FOR  EACH  VIEW  ANGLE  AND  LAYER 
C  PERMUTAION 
C 

PHIT (I,M)»PHIT1 (I,M) 

C 

C.. .CALCULATE  THE  PROBABILITY  OF  GAP  (PGAP)  FOR  ALL  PERMUTATIONS. 

C 

PGAP< I ,M)*1 .-PHIT ( I ,M) 

C 

C . . . CALCULATE  THE  PROBABILITY  OF  GAP  AND  HIT  FOR  THE  HALF  LA YERS ( PGAP2 , PHI T2 ) 
C  FOR  ALL  PERMUTATIONS. 

C 

PGAP2(I,M)=  S0RT(P6AP( I ,M) ) 

PHIT2(I,K)=1.-PGAP2(I,M) 

10  CONTINUE 
C 

C... OBTAIN  THE  FREQUENCY  OF  OCCURENCE  (FREQD)  OF  ELEMENTS  IN  EACH  OF  THE 
C ...  NINE  INCLINATION  INTERVALS  FOR  EACH  LAYER. 

C 

DO  15  J=1 , 3 
ADD=0.0 
DO  20  N=1 ,9 

FREQDtN, J)=  AYLFA(2*N, J) 

ADD=ADD  ♦  FREODIN, J) 


AT  6 


/ 


20  CONTINUE 
DO  25  K*1 ,9 

FREQDIK, J)=FREQD<K,  J)/ADD 
25  CONTINUE 
IS  CONTINUE 
RETURN 
END 

SUBROUTINE  SRVCMQD 
C 
C 

C... SUBROUTINE  SRVCMQD  IS  A  MODIFIED  VERSION  OF  A  PORTION  OF  THE  SRVC 
C  MODEL  WHICH  CAICUALTES  THE  GEOMETRIC  PARAMETERS  OF  A  CANOPY. 

C 

C 

COMMON/StO/  PHI T 1  <3, 9) ,FLAI <3, 1 ) ,SLAI (3, 1 ) ,AXLFA( 19,3) , AYLFA< 1 9,3) 
DIMENSION  N ANGLE ( 3,3) ,FLA<  3,3, 1 0) ,THETA( 10) 

DIMENSION  PHIT(3,3, 10) ,MTP(3) ,OPM( 10) ,XK(9) , XLFA( 1 9) 

DIMENSION  YLFA( 19),  DM(!7),  F(19),  0P(9) 

REAL  INCLF 
r 

C 

C.... GENERAL  SIMULATION  CONSTRAINTS  SRVC 

C 

CEPIQ2=  1 .57079632 
CE2PI=  6.28316530 
CEIPI-  3.1<f592<i5 
CtDTR=. 017453293 
CERTD=  57.2957795 
CEMTR=  .00029088821 
N 6AKBS=9 
NMATs) 

NLAY=3 


BANDU=90/NBANDS 

SRVC 

PARAMETER  INITIALIZATION 

AND  CONVERSION 

SRVC 

NS0UR=NBANDS+1 

SRVC 

BANDW=8ANDW*CEDTR 

SRVC 

COEFFICIENTS  FOR  DIFFUSE 

RADIATION  VECTORS 

SRVC 

ALPKA2=0. 

SRVC 

SINA2=0. 

SRVC 

DO  2  1=1 , NBAKDS 

SRVC 

SINAI =SINA2 

SRVC 

AT  7 


/ 


ALPHA2=ALPHA2+BANDW  SRVC 

SINA2=SIN< ALPHA2)  SRVC 

XK( I )=SINA2*SINA2-S1NA1 *SINA1  SRVC 

2  CONTINUE 
C 

C.... SOURCE  DIRECTION  INCLINATION  ANGLES  SRVC 

C 

TOT AL«0.  SRVC 

THETA ( f )=(BANDW/2. )-BANDU  SRVC 

DO  3  1=1 .NBANDS  SRVC 

THETA(I+i)=THETA(I )+BANDW  SRVC 

3  CONTINUE 
C 

C _ CANOPY  GEOHETRY.  EACH  CANOPY  LAYER  IS  COMPOSED  OF  ONE  OPTICAL  SRVC 

C _ MATERIAL  WHICH  MAY  BE  SPECIFIED  AND  UNIQUE  GEOMETRICAL  PROPERTIES.  SRVC 

C.... CANOPY  GEOMETRIC  PARAMETERS  CONSIST  OF  (l)LEAF  ANGLE  FREQUENCY  SRVC 

C... .DISTRIBUTION  FUNCTION  DENOTED  BY  XLFA  AND  YLFA  ( 2 > LEAF  AREA  INDEX  SRVC 

C _ DENOTED  BY  F LAI  AND  <3)CAN0PY  DENSITY  DENOTED  BY  SLAI.  XLFA  (DEG)  SRVC 

C....AKD  YLFA  MUST  BE  SPECIFIED  AT  AN  ODD  NUMBER  (NANG)  OF  EVENLY  SPACED  SRVC 

C _ POINTS.  FLAI  IS  NON-NEGATIVE  AND  SLAI  RANGES  BETWEEN  0  AND  1.  SRVC 

C 

DELF=10.*CEDTR  SRVC 

DO  350  IL=1 ,NLAY  SRVC 

NANG=1? 

C 

C.. .ASSIGN  THE  NUMBER  OF  MATERIALS  IN  ANY  GIVEN  LAYER 


IhAT=1 

KTP(IL)  =  IMAT  SRVC 

IhATI =IMAT 

DO  351  J= 1 , IMAT 1 

IMAT  =  J 

DO  41  MM=1 ,NANG 

XLFA(MM)=AXLFA(MM,IL> 

YLFA ( MM )=A YLFA (MMfIL) 

41  CONTINUE 
C 


C... .INTEGRATE  AND  NORMALIZE  THE  LEAF  ANGLE  FREQUENCY  DISTRIBUTION  SRVC 

C _ FUNCTION  USING  SIMPSONS  RULE  —  THIS  IS  TEMPORARILY  DENOTED  BY  F.  SRVC 

C....M-1  EQUALLY  SPACED  INTERVALS  OF  F  ARE  THEN  DETERMINED  AND  DENOTED  SRVC 
C ....  BY  FLA  (M  POINTS).  THE  TABLE  FLA  IS  USED  FOR  RANDOMLY  SELECTING  SRVC 

C _ LEAF  INCLINATION  ANGLES.  SRVC 

C 

DO  305  1=1, NANG  SRVC 

305  XLFA(I)=XLFA(I)*CEDTR  SRVC 
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H=< (NANG-1 )/2)+1  SRVC 

NANGLE ( IL, IHAT )=H  SRVC 

CALL  TBLR(M,XLFA,YLFA,BK,F)  SRVC 

DO  310  IANG=1,M  SRVC 

310  FLA(IL,IHAT,IANG)=DH(1ANG>  SRVC 

C 

C 

C. . . .NORMALIZE  THE  INPUT  LEAF  FREQUENCY  DISTRIBUTION  FUNCTION  TO  OBTAIN  SRVC 

C _ A  DENSITY  FUNCTION  F  UHICH  IS  SPECIFIED  AT  ft  POINTS.  SRVC 

C 

FT0T=0.  SRVC 

DO  311  1=1 , NANG  SRVC 

311  FTOT=FTOT+YLFA< I )  SRVC 

DO  312  1=1,9  SRVC 

312  F(I)=(YLFA(2*I)+YLFA(2*I+1))/FT0T  SRVC 

DO  315  1=1, NANG  SRVC 

315  XLFA(I)=XLFA(I)*CERTD  SRVC 

M=M-1  SRVC 

C 

C.... CALCULATE  THE  MEAN  PROJECTION  (OP)  IN  THE  DIRECTION  OF  THE  SOURCE  SRVC 
C... 'THETA)  OF  ONE  UNIT  LEAF  AREA  WITH  INCLINATION  INCLF.  THE  LEAVES  SRVC 

C _ AT  THIS  ANGLE  ARE  ASSUHED  TO  BE  AZIHUTHALLY  ISOTROPIC.  SRVC 

C 

DO  330  IANGLE=1 ,NSOUR  SRVC 

INCLF=-5.*CEDTR  SRVC 

DO  320  1=1,9  SRVC 

INCLF=INCLF+DELF  SRVC 

320  CALL  CQP(INCLF,THETA(IANGLE),0P(I),CEPI02) 

r 

C.... CALCULATE  THE  KEAN  PROJECTION  (OPH)  IN  THE  DIRECTION  OF  THE  SOURCE  SRVC 
C.... (THETA)  OF  ONE  UNIT  LEAF  AREA  AVERAGED  OVER  THE  CANOPY  LEAF  ANGLE  SRVC 

C _ DENSITY  FUNCTION  F.  SRVC 

C 

CALL  COPM(F,OP,OPH(IANGLE) )  SRVC 

C 

C... .CALCULATE  THE  PROBABILITY  OF  A  HIT  (PHIT)  FOR  A  LIGHT  RAY  WITH  SRVC 

C _ SOURCE  DIRECTION  THETA.  SRVC 

C 

CALL  PDENS(IL,IHAT, IANGLE, OPH( I ANGLE), THETA, NANGLE, FLA, SLAI,FLAI, 

*  PHIT) 

330  CONTINUE  SRVC 

351  CONTINUE 

350  CONTINUE  SRVC 

J=NHAT 

DC  228  1=1,3 
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DO  22 8  H*1,9 

PHIT1 (I,M)«PHIT(I,1 ,H+1 ) 

228  CONTINUE 
RETURN 

END  SRVC 

SUBROUTINE  COP (ALPHA, BET A, OP, CEP  1 02) 

C  COP 

C 

C. . . .THIS  PROGRAM  CALCULATES  THE  HEAN  PROJECTION  OF  A  UNIT  LEAF  AREA  IN  COP 

C _ THE  DIRECTION  OF  THE  SOURCE.  THE  LEAF  IS  INCLINED  AT  AN  ANGLE  COP 

C... .ALPHA  AND  IS  ASSUMED  TO  BE  AZIMUTHALLY  ISOTROPIC. 

C  COP 

C 

C 

OPrCOS (ALPHA) *SIN( BETA)  COP 

IF(ALPHA.LE.BETA)  RETURN  COP 

C 

C....THETAO  IS  THE  LEAF  AZIMUTH  ANGLE  AT  WHICH  OP  BECOMES  NEGATIVE  AND  COP 

C _ IS  IN  THE  FIRST  QUADRANT.  THE  FUNCTION  OP  IS  SYMMETRIC  AND  HENCE  COP 

C _ IS  AVERAGED  OVER  LEAF  AZIMUTH  ANGLES  OF  0  TO  PI  RADIANS.  COP 

C 

THETAO=ACOS ( TAN (BETA) /TAN( ALPHA))  COP 

TANTO=TAN(THETAQ)  COP 

OP=OP*< I .+(TANT0-THETA0)/CEPIC2)  COP 

RETURN  COF 

END  COP 

SUBROUTINE  COPM(G,OP,OPM)  CORN 

C 
C 

C. . . .THIS  PROGRAM  CALCULATES  THE  MEAN  PROJECTION  OF  A  UNIT  LEAF  AREA  IN  COPM 
C....THE  DIRECTION  OF  THE  SOURCE  (OPM)  FOR  THE  SIMULATED  CANOPY.  THE  COPh 

C _ LEAVES  OF  THE  CANOPY  ARE  ASSUMED  TO  BE  AZIMUTHALLY  ISOTROPIC.  THE  COPK 

C....OP  FUNCTION  USED  IN  THE  CALCULATION  HAS  BEEN  PREVIOUSLY  DETERMINED  COPM 

C. . . .FOR  A  GIVEN  SOURCE  DIRECTION  FOR  LEAF  INCLINATION  ANGLES  OF  COPM 

C _ 5,  13,  ...,  85  DEGREES.  G  IS  THE  LEAF  INCLINATION  ANGLE  DENSITY  COPM 

C... .FUNCTION.  COPM 

C  COPM 

C 

DIMENSION  0P(9) ,G(9)  COPM 

0PM=0.  COPM 

DO  1  1*1,?  COPM 

1  OPMsOPM+OP( I )*G( I )  COPM 

RETURN  COPM 

END  COPM 

SUBROUTINE  PDENS< IL,HTYPE, I  ANGLE, OPM, THETA, NANGLE, FLA  ,SLAI  ,FLAI , 
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*  ph:t> 


c 

c 

C - THIS  PROGRAM  COMPUTES  THE  PROBABILITY  THAT  LIGHT  AT  INCIDENT  ANGLE  PDENS 

C  T  HE  T  A  < I  ANGLE )  INTERACTS  WITH  MATERIAL  TYPE  MTYPE  U1THIN  CANOPY  PDENS 

C  LAYER  IL.  PDENS 

C  PDENS 

C 

C  INPUT  PDENS 

C  IL  PDENS 

C  MTYPE  PDENS 

C  IANGLE  PDENS 

C  OPM  PDENS 

C  SLAI  PDENS 

C  FLA I  PDENS 

C  THETA  PDENS 

C  OUTPUT  PDENS 

C  PHIT  PDENS 

C  PDENS 

C 

DIMENSION  DUM<357),THETA(10) 

DIMENSION  NANGLE<3,3),FLA<3,3,10),SLAI<3,3>,FLAI(3,3),PHIT(3r3,10)  PDENS 
ARG=1 . -(SLAI ( I L, MTYPE )*OPM/SIN<  THETA ( I  ANGLE ) ) )  PDENS 

IF  (ARG.LE.O.)  GO  TO  1  PDENS 

PO*ARG**(FLAI ( IL,MTYPE)/SLAI < IL f MTYPE ) )  PDENS 

GO  TO  2  PDENS 

1  PO  *  0.  PDENS 

2  CONTINUE  PDENS 

PHITtIL, MTYPE, IANGLE)=1 ,-PO  PDENS 

RETURN 

END  PDENS 

SUBROUTINE  TBLR<M,  X,  Y,  XX,  Z)  TBLR 

C  TBLR 

C 

C....THIS  PROGRAM  FINDS  THE  INTEGRAL  Z(X>  OF  THE  FUNCTION  Y ( X )  FROM  X ( 1 )  TBLR 

C _ TO  X(2H-1>  USING  SIMPSONS  RULE.  THE  INTEGRAL  Z(X)  IS  NORMALIZED  TO  TBLR 

C _ 1.0  AT  X ( 2H- 1 ) .  THE  TABLE  OF  Z  VERSUS  X  IS  THEN  INVERTED  TO  DETER-  TBLR 

C _ MINE  X  AS  A  FUNCTION  OF  Z  AT  M  REGULARLY  SPACED  POINTS  ALONG  Z.  TBLR 

C  TBLR 

C 

C  INPUT  VARIABLES  TBLR 

C  «  =  DESIRED  NUMBER  OF  REGULARLY  SPACED  POINTS  ALONG  Z  TBLR 

C  X  =  SPECIFIED  AT  2M-1  POINTS  TBLR 

C  Y  =  SPECIFIED  AT  2M-1  POINTS  TBLR 

C  OUTPUT  VARIABLES  TBLR 
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t 


I 


i 

i 

|! 

f 


c 

XX  *  THE  TABLE  OF  X  VALUES  FOR  M  REGULARLY 

SPACED  POINTS 

TBLR 

c 

(M-1  INTERVALS)  ALONG  Z. 

TBLR 

c 

Z  *  THE  NORMALIZED  INTEGRAL  OF  Y  AT  X(1>, 

X<3),  ...,  X(2M-1 ). 

TBLR 

c 

c 

c 

c... 

r 

TBLR 

DIMENSION  X(19),Y(I9),2(10),XI(10),XX(10) 

TBLR 

.SIMPSONS  RULE  INTEGRATION 

TBLR 

l 

10 

Z<1)  *  0.0 

TBLR  }: 

DX  *  X  <  2 )  -  X  <  1  > 

TBLR 

20 

DO  50  J  «  2,M 

TBLR 

JO  =  2*J  -  3 

TBLR 

30 

J1  *  2*J  -  2 

TBLR  i 

J2  *  2*J  -  1 

TBLR  ; 

40 

Z(J)  =  Z(J  -  1)  +  DX<( Y( JO)  +4.*Y(J1)  +  Y(J2)  )/3.0 

TBLR 

50 

XI(J>  *  X ( J2 ) 

TBLR  { 

XI(1)«X<1) 

TBLR  | 

w 

C  •  •  • 

r 

.NORMALIZE  INTEGRAL  Z(X) 

TBLR  * 

L 

60 

DO  70  J  =  l ,M 

TBLR 

70 

r 

Z(J)  =  Z( J)/Z(M) 

TBLR 

t 

C  •  •  • 

r 

.FIND  X  AT  M  REGULARLY  SPACED  POINTS  ALONG  Z. 

TBLR 

i 

V 

X 

11 

X 

X 

TBLR 

EM  *  M  -  1 

TBLR 

F  *  1.0/EM 

TBLR 

JS=2 

TBLR 

80 

DO  120  K  =  2,N 

TBLR 

ZT  =  K  -  1 

TBLR 

ZT  «  ZT*F 

TBLR 

?0 

DO  110  J  =  JS,H 

TBLR 

I F < Z ( J )  -  ZT)  110,  100,  100 

TBLR 

100 

G  =  (ZT  -  Z(J  -  1>)  / tZ( J>  -  Z(J  -  D) 

TBLR 

XX(K)  «=  XI  (J  -  1)  +  G*(XI(J)  -  XI  ( J  -  D) 

TBLR 

GO  TO  115 

TBLR 

110 

CONTINUE 

TBLR 

115 

JS  =  J 

TBLR 

120 

CONTINUE 

TBLR 

RETURN 

TBLR 

END 

SUBROUTINE  SECTOR 

TBLR  i 

i 

c 


s 

i 

i 


c 

C... SUBROUTINE  SECTOR  CALCULATES  THE  NORMALIZING  FACTORS  WHICH  ACCOUNT  FOR  THE  A 
C  AREA  OF  EACH  SOURCE  SECTOR. 

C 

C 

COMMON  /A/  UV,RH,RL,D<3) 

COMMON  III  P6AP<3,9>,  PHIT<3, 9 ) ,PGAP2(3,9 > ,  PHI T2 ( 3 r 9) ,STEF 
COMMOK/C/COSTA( 9 ,9 , 1 6) , EMISSVY 4 ) , ABSORB ( 3 > ,  ESKY,  EGRD,  SECTARY) 
COHMON/D/COKT (3,5,9) ,C<3,5,9) ,  SUMT(3,9),  KELV,  GT,  NUSIM,  ITIME 
COMMON/E/ AT ,THETA( 9) , PHI ( 1 8) ,XLF(9) ,  YLF ( 9) ,  ZLF(9),  XS(9,t8) 
C0MM0N/F/YS(9,16),ZS(9),  CEDTR,  B,  FRECD(9,3>  ,  UA(15),EPS 
C0MM0N/G/NSIG,N,  UMAX 
BANDW*  IO.*CEDTR 
ALPHA2®  0. 

SINA2*0 . 

DO  2  1=1,9 
SINAI *SINA2 

ALPHA2®  ALPHA2  ♦  BANDU 
SINA2=  SIN  (ALPHA2) 

C 

C...  NDTE  UE  MUST  DIVIDE  BY  SIN(THETA)  SINCE  UE  ARE  INTERESTED  IN  THE  FLUX 
C  BEFORE  IT  HITS  A  HORIZONTAL  PANAL. 

C 

SECTARY  I >=  (SINA2»*2-SINA1**2)/(18.*SIN(THETA(I))) 

2  CONTINUE 
RETURN 
END 

SUBROUTINE  SET02Y A,  I , J) 

C 

C 

C.. .SUBROUTINE  SET02  SETS  ALL  ELEMENTS  OF  A  2-DIMENSIONAL  ARRAY  TO  0.0 

C 

C 

DIMENSION  A< I , J ) 

DO  10  K«1,I 
DO  10  L=1 , J 
A < K , L ) *  0.0 
10  CONTINUE 
RETURN 
END 

SUBROUTINE  SET03  (A,I,J,K) 

C 

C 

C... SUBROUTINE  SET03  SETS  ALL  ELEMENTS  OF  A  3-DIMENSIONAL  ARRAY  TO  0.0 
C 
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DIMENSION  AU'J'K) 

DO  tO  L*t , I 
DO  tO  M*1,J 
DO  tO  N«1 ,K 
A(L,M,N)«  0.0 
tO  CONTINUE 
RETURN 
END 

SUBROUTINE  UMAT (TITLE ) 

CONNON/B/  PGAP<3,9),  PHI T < 3 , 9 ) ,  PGAP2<3,9>,  PHIT2<3,9) , STEF 
DIMENSION  U(4,9),TITL£(8> 

C  UNAT  CALCULATES  THE  U  MATRIX. 

C 

DO  10  M  =  1,9 
U( 1  ,M)  >  PHIT< 1 ,M) 

U(2,M)  =  PGAP ( 1 ,M)-P6AP< t,M)*PGAP(2,M) 

U(3.K)  =  PGAP( 1 ,M)*PGAP<2,M)-PGAP< 1 ,M)*PGAP(2,M)*PGAP<3,H) 
UU,M>  *  PGAP( 1 ,M)*PGAP(2,M)*PGAP(3,N) 
tO  CONTINUE 

URITE(6, 1 99) (TITLE(N) ,N*1 , 8) 

URITE(3,199)(TITLE<N),N=t,8) 

199  FORMAT (*  THE  U  AND  S  MATRICES  FOR  */8A10) 
URITE(3,300)(<W<M,N),N=t ,9),M=1 ,4) 
URITE(6,300)<<WIM,N),N=1,9),H=1,4) 

300  FORMAT (9F7.4 ) 

RETURN 

END 
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SENSIT 


PROGRAM  SENSIT ( INPUT, OUTPUT, TAPES* INPUT ,TAPE6a0UTPUT,TAPE1 ,TAPE2) 

CQNN0N/PARA1/S1Q,STAC3),S<3,3>,STGC3),X<3),PARAN<1A),AC3>, 
1ALP(3),FNU,TA,TAC,T6,EPSC3),EPSTG,RH,RL 
CONNQN/PARA2/BTA,BTG,3X( 3) , DBX ( 3) , 8XC  3) , DQXC  3) ,RX(3) , DRX ( 3) 
CQHN0N/ESTIM/FXC3) , DFX(3,3) 

DIMENSION  PAR(  14),  DSTEPC 1 6) ,DP(  1  A) ,DX(3) , TENPC3) ,X0(3) , XXC  3) , 
tFSEMC3>,PARMAXCU> 

— -  DEFINE  INPUT  VARIABLES  — 

SI3a5. A484E-8 
STAC1)a0,4107 
STA(2)*0. 1 887 
STAC3)*0.0728 
STG(t>*0.0815 
STGC2)a0.2482 
STG(3)=0.7257 

SCI  ,1  ) *0.7715  $  S(1,2)a0.47A9  »  SC  1  ,3)=0.0523 

S(2,1  )*0.22 77  i  S(2,2)=l .1600  i  S<  2, 3)»0. 1 682 

SC3,1)a0.0830  *  SC  3,2 )a0. 5498  «  SC3,3)=0.5414 

—  DEFINE  NOMINAL  VALUES  FOR  MODEL  PARAMETERS  - 

PAR( 1 )aALP< 1 )*1 .0 
PAR(2)=ALPC2)a1.Q 
PAR<3)»ALPC3>a1.0 
PAR(4)»FMU»10.0 
PAR(5)»TA«14.6 
PAR<6)»TACa14.4 
PAR<7)*TG*1 1 .7 
PARC  8)»EPSC 1 )*1 .0 
PARC9)»£PSC2>a1.0 
PAR(10)aEPS(3)a1 .0 
PARC  1 1 )*EPSTG*1 .0 
PAR(t2)aRH»0.20 
PARC  1 3)*RLa0. 44 
PARC  1 4)*A< I )a1 44. 

PARn5)»AC2)«49, 

PARC  1 A)aAC3)a4A. 

—  DEFINE  INITIAL  VALUES  FOR  STATE  VARIABLES  — 

X0< 1 )»20.4 
X0(2>a14.A 


n  o  o 


X0<3)*16.3 


—  DEFINE  UPPER  BOUNDS  FOR  NODEL  PARAMETERS  — 

PARHAX< 1 )»PARNAX(2)*PARMAX(3)*1 .0 
PARNAX! 8) -PARNAX! 9) -PARNAX! 1 0)=PARNAX ( 1 1 )  *1 .0 
PARNAX(4)*20.0 

PARNAX !5)-PARNAX!6)*PARNAX(7)-30.0 
PARNAX! 12)-0. 50 
PAR«AX<13>»1.20 
PARNAX ( 1 4)-1 50. 0 
PARNAX! 15J-PARNAX! 16)*60. 

C 

C  —  DEFINE  STEP-SIZES  OF  NODEL  PARAMETERS  FOR  SENSITIVITY  ANALYSIS 
C 

DSTEPt 1 )*DSTEP<2)=DSTEP(3)=DST£P!3>-DSTEP(9)=DSTE?< 1 0)*DSTEP(  1 1 ) » 
(-0.005 
DSTEP<4)=- 1.0 

DSTEP!5)»DSTEP<6)=DSTEP<?>*-1 .0 
DSTEPU2)— 0.05 
DSTEPt 1 3>»- 0.05 
DSTEP!  1 4)— 7.5 
DSTEPt I5)*DSTEP< 1 A)»-3.0 
C 

C  —  START  SENSITIVITY  ANALYSIS  - 

C 

TOl-O.OOOl 
DO  100  IP-1,16 
DO  92  1-1,16 
92  PARAN!I>«PAR!I) 

DO  94  1-1,3 

xxm-xom 

94  XU)-XO(I) 

DO  90  NTINE-1 ,20 

FACTOR-NTINE-1 

DP! IP)-FACTOR*DSTEP( IP) 

PARAN!IP>«PARNAX<IP)+DP!IP> 

50  CALL  FEVAL 

C  WRITE(6,1 010MFX!  I  >  ,1-1 , 3) 

C1010  FORMAT! 1H  ,3F10.5) 

DO  20  1-1,3 
20  F X <  I )— FX( I ) 

CALL  SQLVE<FX,3,DFX,3,DX) 

DO  30  1-1,3 

C  WRITE ( 6, 1  01 0)  X(I),DX(I) 
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30  X<I)aX<I)+DXU> 

DO  40  1*1,3 
D£V=DX( I ) 

IF < A3S( DEV)  .ST.  TOL)  SO  TO  50 
F3EM(I)*<X(I)-XX(I))/DSTEP(IP) 

XX(I)»X(I) 

40  TEMPI  I )  =  X  1 1 ) 

URITEM ,90)  ( TEMPI  I ) ,1*1 ,3) ,?ARAH< IP) 

30  FORMAT ( 3F 1 0 .5) 

90  CONTINUE 
100  CONTINUE 
SUBFILE  1 
REUIND  1 
STOP 
END 

SUBROUTINE  FEUAL 

COMMON/PARA1 /SIB,STA(3) ,S<3,3),ST3(3),X(3),PARAH(16),A(3>, 

1  AL?<  3) ,FMU, TA,TAC,TG,EPS<  3) ,EPSTG,RH,RL 
COMMON/ PAR A2/3TA,3TG,3X (3) ,DBX( 3) ,9X(3) ,DQX(3) ,RX(3) ,DRX(3) 
C0MMQN/£STIH/FX(3),BFX(3,3> 

ALP( 1 }*PARAN< 1 ) 

ALP(2)=PARAN(2) 

AL?(3)*PARAM<3) 

FJ1U=PASAM(4) 

TA»PARArt(5) 

TAC=PARAM<6) 

TG*PARAN(7) 

EPS11 )=PARAM<8) 

EPS(2)*PARAM<9) 

EPS(3)=PARAX<10) 

E?STG*PARAM(  1 1 ) 

RH=PARAM< 12) 

RL*PARAM( 13) 

A< 1 )=PARAM< 1 4) 

A(2)*PARA«( 15) 

A(3)*PARAM( 14) 

CALL  F3TA(TA,3TA) 

CALL  DFUNC(EPSTG,TG,DTG,D3TG) 

DO  10  1=1,3 

CALL  3FUNC(EPS( I) ,X( I) ,3X( I ) , D3X( I ) ) 

CALL  OFUNC(X(I),TAC,F»U,OX(I),DaX(I)) 

CALL  RFUNC (X( I )  ,TAC,FMU,RL,RH,SX< I ) , DRX( I ) ) 

10  CONTINUE 
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DO  20  IL*1 ,3 

20  FX<IL)=0.5*AL?<lL>*SIG*<BTA*STA<IL)+3X(1>*saL,1)+BX(2>* 
1S<IL,2>+BX(3)*SCIL,3)  +  3TG*S7G<im+A<IU-SIG*3X(IL>+aX<Il> 
2+RX(IL) 

DO  30  1*1,3 
DO  30  J*1 ,3 
30  DFX< I , J)=0. 

DO  40  1L*1 ,3 
DO  40  J*t  ,3 
IFU.NE.IL)  GO  TO  35 

DFX(IL,J>=0.5*ALP<IL)*SI3«D2X(J)*$(IL,J>-3IG*D3X(IL)+DQX<IL>+ 
1 DRX< IL) 

GO  TO  40 

35  DFX1IL, J)*0.5*ALP(IL)*SIS*D8X( J)*S(IL, J) 

40  CONTINUE 
RETURN 
END 

SUBROUTINE  BF UNC < EPS I , X I , BX I , D3X I ) 

BXI*£P$I*<XI+273.0)**4. 

DBXI*4. *EPSI*( XI +273.0 )**3, 

RETURN 

END 

SUBROUTINE  F3TAdA,3TA) 

EPSTA*1 .-0.261 *£XP<-?.77E-4*TA*TA) 

CALL  BFUNC(EPSTA,TA, BTA, BBTA) 

RETURN 

END 

SUBROUTINE  OFUNC(XI , TACfFNU,QXI ,IiQXI ) 

IF(FNU.GT.30. )  GO  TO  10 
HC=0.69775*(20.4+0.2*FHU**.97) 

GO  TO  20 

10  HC=0. 69775* (0.95*FHU**. 97) 

20  QXI=(XI-TAC)*HC*(-1 .0) 

D0XI=HC*(-1 .0) 

RETURN 

END 

SUBROUTINE  RFUNC(XI,TAC,FHUrRL,RH,RXIf DRXI) 

RNUM=FEX(XI)*1.0E-6-RH*FEX(TAC)*1 .OE-6 
RDEN=RL  +  ( 1 ./60. )* <0.04  +  1  .27*F«U**(-0.5) ) 

RXI1 =-697.75*(-0.566*XI+597.3) 

RXI2=RNUH/RDEN 

RXI*RXI1*RXI2 

DRXI=697. 75*0. 566*RXI2+RXI1*( 0.05671 5E-6*FEX( XI ) )/RBEN 

RETURN 

END 
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SUBROUTINE  1NVERSE1 A,M,D) 

—  INVERT  A  3*3  REAL  NATRIX  A  WHOSE  DETERMINANT  IS  D 
THE  RESULT  U1LL  BE  STORED  IN  A 

DIMENSION  A(3,3),B<3,3> 

D*A<1,1)*A<2,2>*A<3,3)+A<1,2>*A<2,3>*A(3,1>+A<1,3)*A<2,1)* 
1A(3,2)-A(3,1)*A(2,2)*A<1,3)-A(I,1)*A(3,2)*A(2,3)-A(2,1)* 
2A( I ,2)*A(3,3) 

B(1,1)a(A(2,2)*A(3,3)-A(2,3)*A(3,2))/D 
B(1,2)*-<A(2r1)*A(3f3)-A<2f3)*A<3,1))/D 
B(1,3)*(A(2,t)*A(3,2)~A<2,2)*A(3,1))/D 
B<2, 1 )*-< A( 1 ,2)*A<3,3)-A( 1 ,3>*A(3,2) )/D 
B(2,2)»(A(1f1)«A(3,3)-A(1,3)*A(3f1))/D 
B(2,3)»-(A<1,1)*A<3,2>-A<1,2)*A<3,1))/D 
B(3,1)*(A(1,2)*A<2,3)-A<1,3>*A(2,2)>/D 
B(3,2)*-<A(1 ,1 )*A<2,3)-A(1 ,3)*A(2,1 ))/D 
B(3,3)*(A(1 ,1 >*A<2,2)-A(1f2)*A(2,1  ))/D 
DO  10  1*1 ,3 
DO  10  J*1 , 3 
10  AU,J)»B(I,J> 

RETURN 

END 

FUNCTION  FEX < XI ) 

XX»3.2342*EXP(0.054?15*XI> 

FEX*XX 

RETURN 

END 

SUBROUTINE  SOLVE! Y, N, A,«,X) 

DIMENSION  Y<N)fA(N,M),X(M>,ATA<M,M) r ATY(M) 

DIMENSION  Y(3) ,A(3,3),X(3),ATA(3,3),ATY<3) 

DO  10  1*1, N 
DO  10  J*1,M 
ATA( I, J)*0. 

DO  10  K*1 f N 

10  ATA( I,J)*ATA<I, J)*A<K,I)*A(K, J) 

CALL  INVERSE! ATA,N, D> 

DO  20  1*1, M 
ATY(I)*0. 

DO  20  J*1 ,N 

20  AT Y ( I ) *AT Y (I)+A(J,I)*Y<J) 
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DO  30  1*1, n 

xm«o. 

DO  30  J*1,N 

30  X<I>*X<I>+ATA<It J>*ATY<J> 
RETURN 
END 


♦DECK  SRVC 

PROGRAM  SRVC< INPUT, OUTPUT, TAPE6»0UTPUT, TAPES" INPUT) 

C....  SOLAR  RADIATION  -  VEGETATION  CANOPY  REFLECTANCE  MODEL 

C....  THIS  PROGRAM  CALCULATES  THE  APPARENT  DIRECTIONAL  REFLECTANCE  OF  A 

C....  VEGETATION  CANOPY  AS  A  FUNCTION  OF  CANOPY  GEOMETRY,  LEAF  REFLEC- 


c. 

...  TANCE 

AND  TRANSMISSION.  SOIL  REFLECTANCE,  AND  CANOPY  IRRADIAMCE 

SRVC 

c. 

...  FOR  A 

GIVEN  SOLAR  POSITION. 

SftVC 

c. 

•  • •  R.E . 

OLIVER  AND  J.A.  SMITH  COLORADO  STATE  UNIVERSITY  JUNE, 

1974 

SRVC 

c 

SRVC 

c. 

SRVC 

c 

SRVC 

c 

LABEL 

EXTERNAL  REFERENCES 

SRVC 

c 

SRVC 

c 

Cl 

BLOCK  DATA,  LANBTN,  SUN,  ETHRES,  LANGLE,  NRH,  SETZ, 

UTIL, 

SRVC 

c 

AND  COP. 

SRVC 

c 

SRVC 

c 

C2 

LANBTN,  PDENS,  AND  OPTICAL. 

SRVC 

c 

SRVC 

c 

C4 

LAN6LE,  PDENS,  AND  PGAP. 

SRVC 

c 

SRVC 

c 

Ci 

ETHRES,  SETZ,  AND  LANBTN. 

SRVC 

c 

SRVC 

c 

C3 

LANGLE. 

SRVC 

c 

SRVC 

c 

Lt 

OPTICAL. 

SRVC 

c 

SRVC 

c 

CHAT 

PGAP  AND  LANBTN 

SRVC 

c 

SRVC 

c 

// 

LANGLE  AND  LANBTN. 

SRVC 

COMMON/C1 /DAY, YEAR, TIME, GLAT,GL0N6, DEC, BANDU,NLAH,THETS1 ,THETS2, 

1 NHAT,EXTRA(4) ,NQP, INIT,DUM1 (13), 

2CEDTR,CERTD,CEMTR,CEPI02,CE1PI ,CE2PI ,BUM2(  1  A) , 

3SINLAT, COSLAT, SINDEC,COSDEC, COSH, SINZ,COSZ,$INAZ,COSAZ,LXS,LYS,LZS 
C0MM0N/C2/CANRMU7)  ,SKYIM<  17), DIFIH(17),R(17),TU7),RG  (17), XL  AM(17 
t ), SOURCE <  10,17) ,THETA( 1 0) ,ZENITH< 10) 

C0MM0M/C4/NANGLE(3,3) ,FLA<  3,3, 1 0) ,SLAI (3,3) ,FLAI <3,3>,PHIT<3,3,10) 
COHHON/CA/DR(4, 10, 17), UR (A, 10, 17), THRESD< 10), I  SOD  <  4,10) , I SOU ( 4, 10) 
1 ,  THRESU  < 1 0 ) 

COMMON  /KIN/  I ML <3,3,2) 

C0NH0N/C3/SINL,C0SL,SINP,C0SP 

C0MMQN/L1/DATAID<7,9) ,XNU( 17,?),C(17,17,?) ,NVEC(9) 
C0HH0N/CNAT/HTP(3),NLAY,0PH<  10) 

A,ENDLC 

COMMON  AV£C(17),XX(9),SXL,SYL,SZL,XLF,YLF,ZLF 
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1,XSU0,19>,YS(10,18>,ZS<tO) 

A,ENDBB 

COMMON  /AB3/TABS0<4,17> 

C . .  INTERNAL  ARRAYS  . 

DIMENSION  JOB  I D  <  8  >  r  VEC  T  <17),S1G<)7),V(17,I7) ,  COR < 1 7, 1 7) 
DIMENSION  COV<  1 0, 17, 17), C0VNM7, 17) 

DIMENSION  XIFA< 1 9) , YLFA( 1 9) ,DN< 1 7) ,DHt (17), REFER (17) 
DIMENSION  RIT (10,17) r R I T BAR ( 1 0, 17) ,RBAR( 1 0, 17) 

DIMENSION  F(19),0P<9) 

INTEGER  RORT 
REAL  LXS,LYS,LZS,INCLF 
INTEGER  DAY, YEAR, TH,TM,ZDEG 
8000  CONTINUE 

DO  10  1*1,10 
THETA  < I >  *0 . 

ZENI7H(I)*0. 

ZS(I)*0. 

THRESD ( I ) *0 . 

0PM(I)*0. 

THRESU(I)*0. 

DO  10  J»1 ,18 
10  YS<1, J)*0. 

DO  4  K* 1 ,17 
CANRH<K)*0. 

SKYIM<K)*0. 

DIFIM(K)*0. 

R(K)*0. 

T(K)*0. 

RG(K)*0. 

DM(K)*0. 

DHt  <K)»0. 

SIG(K)*0. 

XLAH(K)*0. 

DO  A  1*1,10 
S0URCE1 I,K)*0. 

RIT< I ,K)*0. 

RITBAR(I,X)*0. 

RBAR<I,K)*0. 

A  CONTINUE 
DO  9  1*1,19 

F(I)*0. 

XLFA(I)*0. 

9  YLFA(I)*0. 

DO  12  1*1,9 
NVEC<1)*0. 
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0P< I)-0. 

XK( I )*Q. 

DO  12  J-1 ,17 
12  XMU< J,I)«0. 

DO  7  1*1,3 

NTP<I)*0. 

DO  7  J*1 ,3 
NAN6LE<I,J)*0. 
SLAI< I , J)*0. 
FLAI (I, J)*0. 

DO  7  K»1 ,10 
FLA< I , J,K)*0. 

7  PHIT < I, J,X)*0. 
DO  6  1*1,3 

DO  6  J»1,3 
DO  i  K«1 ,2 
6  INL< I , J,K)*0. 


C... .PERIPHERAL  CONTROLS  SRVC 

I HI ST  *  0  SR VC 

ISTOH  *  1  SR VC 

IFILE  •  5  SR VC 

C... .IFILE  ASSISNNENT  COULD  BE  HADE  THRU  A  READ  STATEMENT. 

IF<EOF<5).NE.O.)  STOP  SRVC 

IF < IHIST.EQ. 1 )  CALL  FUN<-1,-1>  SRVC 

C.... GENERAL  SIMULATION  CONSTRAINTS  SRVC 

R£AD< IFILE, 100)  JQBID,DAY, YEAR,TH,TM,GLAT ,i»LQNG,DEC,NBANDS,  SRVC 

1NLAft,MMA7,IMZT,NSAMP,NTRlAL  SRVC 

IF(EOF<5).NE.O.)  STOP  SRVC 

READ<IFILE,102>  NLAY  SRVC 

BANDU*90/NBANBS  SRVC 

URITE(6,200>  JOBID,DAY,  YEAR,TH,TM,GLAT,0LONG,DEC,BANDU,NLAM,NMAT ,  SRVC 
1 INIT ,NSAMP,NTRIAL,NLAY  SRVC 

READ< IFILE,  1 01 )  THRESD  IREADdFILE,  101 )  THRESU  SRVC 

UR1TE<4,221 )  THRESD, THRESU  SRVC 

C.... PARAMETER  INITIALIZATION  AND  CONVERSION  SRVC 

DO  1073  J*1,4 
DO  1073  1*1 ,NLAM 
TABSO< J,I)*0.0 
1073  CONTINUE 

NSQUR*NBANDS+1  SRVC 

HLAYP1*NLAY*1  SRVC 

CALL  RANSET(IMIT)  SRVC 

XT1*TH  SRVC 

XT2*TN  SRVC 

TIME*XT1  MXT2/A0. )  SRVC 
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GLAT*6LAT*CEDTR  SRVC 

GL0M6*6L0N6*CEDTR  SRVC 

DEC*DEC*CEDTR  SRVC 

BANDU*8ANDW*CEDTR  SRVC 

C. . . .SUN  POSITION  PARAMETERS  SRVC 

CAU  SUN  SRVC 

URITE<6,222)  LXS,LYS,LZS  SRVC 

ZS<  1 )  *  LZS  SRVC 

C.... COEFFICIENTS  FOR  DIFFUSE  RADIATION  VECTORS  SRVC 

C....SENSOR/9AND  AREA  RATIO  FOR  ALL  DIFFUSE  BANDS  SRVC 

ALPHA2*0.  SRVC 

SINA2»0.  SRVC 

DO  2  I*1,NBANDS  SRVC 

SINAI *SINA2  SRVC 

ALPHA2*ALPHA2+BANDU  SRVC 

SINA2*SIN<  ALPHA2)  SRVC 

XX(I)»SINA2*SINA2-5INAI»SINA1  SRVC 

2  CONTINUE 

URITE <6,208)  (XK( I ) ,1*1 ,N8ANDS>  SRVC 

C... .SOURCE  DIRECTION  INCLINATION  ANCLES  SRVC 

TOTALS.  SRVC 

THETA< 1 )*<3ANDU/2. )-BAMDU  SRVC 

DO  3  1*1 , NBANDS  SRVC 

THETA < 1*1 )*THETA( I )+BANDW  SRVC 

3  CONTINUE 

THETA( 1 )*CEPI02-AC0S(CQSZ)  SRVC 

CONS»LZS*TOTAL  SRVC 

DO  50  1*1,10  SRVC 

50  ZENITH( I )*CEPI02“THETA( I )  SRVC 

URITE(6,223)  THETA  SRVC 

C.... DIRECTION  COSINES  OF  AZIMUTHAL  SECTORS  IN  THE  DIFFUSE  BANDS  SRVC 

DEG20*20.*CEDTR  SRVC 

DO  60  JS0R*2,NS0UR  SRVC 

ZS( JSOR) *SIN( THETA (JSOR ) )  SRVC 

PHI»10.*CEDTR  SRVC 

DO  60  IPHI*1 ,18  SRVC 

XS( JSOR, IPMI)*CQS( THETA (JSOR) )*COS(PHI )  SRVC 

YS( JSOR, 1PH1)*C0S(THETA4 JSOR) )*SIN(PHI)  SRVC 

60  PHI«PHI*D£620  SRVC 

C.... CANOPY  SEONETRY.  EACH  CANOPY  LAYER  IS  COMPOSED  OF  ONE  OPTICAL  SRVC 

C.... MATERIAL  WHICH  MAY  BE  SPECIFIED  AND  UNIQUE  GEOMETRICAL  PROPERTIES.  SRVC 

C... .CANOPY  GEOMETRIC  PARAMETERS  CONSIST  OF  (DLEAF  ANGLE  FREQUENCY  SRVC 

C.... DISTRIBUTION  FUNCTION  DENOTED  BY  XLF A  AND  YLFA  (2)LEAF  AREA  INDEX  SRVC 

C... .DENOTED  BY  FLAI  AND  (31CAN0PY  DENSITY  DENOTED  BY  SLAI.  XLFA  (DEG)  SRVC 


C. . . .AND  YLFA  MUST  BE  SPECIFIED  AT  AN  ODD  NUMBER  (NANG)  OF  EVENLY  SPACED  SRVC 
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C... .POINTS.  FLAI  IS  NON-NEGATIVE  AMD  SLAI  RANGES  BETUEEN  0  AND  1 .  SRVC 

DELF-1 0. *CEDTR  SRVC 

UR I T  E  <  6 , 2  27 )  SRVC 

DO  350  IL»t , MLAY  SRVC 

R£ADdFILE,102)  NANG  SRVC 

C. . .READ  IN  THE  NUMBER  OF  MATERIALS  IN  ANY  GIVEN  LAYER 

READdFILE,  1 02)  INAT  SRVC 

MTP < IL >  «  INAT  SRVC 

I MAT  1*1 MAT 
DO  351  J*1,INAT1 
INAT  *  J 

READ( IFILE, 1 01 )  (XLFA< I )  ,YLFAd )  ,1*1 ,NANG)  SRVC 

READdFILE,  101 )  SLAKIL, INAT), FLAKIL, INAT)  SRVC 

C.... INTEGRATE  AND  NORNALIZE  THE  LEAF  ANGLE  FREQUENCY  DISTRIBUTION  SRVC 

C... .FUNCTION  USING  SIMPSONS  RULE— THIS  IS  TEMPORARILY  DENOTED  BY  F.  SRVC 

C....N-1  EQUALLY  SPACED  INTERVALS  OF  F  ARE  THEN  DETERMINED  AND  DENOTED  SRVC 

C....3Y  FLA  <N  POINTS).  THE  TABLE  FLA  IS  USED  FOR  RANDOMLY  SELECTING  SRVC 

C _ LEAF  INCLINATION  ANGLES.  SRVC 

DO  305  1*1, NANG  SRVC 

305  XLFA(1)*XLFA( I )*CEDTR  SRVC 

M*( (NANG-1 )/2) ♦!  SRVC 

NANGLE<IL,IHAT)»M  SRVC 

CALL  TBLR(H,XLFA, YLFA,DM,F)  SRVC 

URITEU,233)  (F(I),I*1,M)  SRVC 

DO  310  IANG*1,M  SRVC 

310  FLA< IL,1MAT,IANG)*DM<IANG)  SRVC 

C.... NORMALIZE  THE  INPUT  LEAF  FREQUENCY  DISTRIBUTION  FUNCTION  TO  OBTAIN  SRVC 
C....A  DENSITY  FUNCTION  F  UHICH  IS  SPECIFIED  AT  M  POINTS.  SRVC 

FT0T*0.  SRVC 

DO  311  1*1, NANG  SRVC 

311  FTOT*FTOT  *YLFA< I )  SRVC 

DO  312  1*1,9  SRVC 

312  F(I)*(YLFA(2*I)*YLFA<2*K1 >)/FTQT  SRVC 

DO  315  1*1, NANG  SRVC 

315  XLFAd  )*XLFA(I  )*CERTD  SRVC 

URITEU,230)  IL,IMAT,  NANG,  <XLFA(I ) ,YLFA< I > ,  1*1,  NANG)  SRVC 

URITE(4,23t )  MANGLE < IL, IMAT )  SRVC 

URITEU,232)  (FLA<IL,IMAT,I),I*1 ,M)  SRVC 

M*M-1  SRVC 

UR I TE < 6 y 233 >  (F(I),I«1,M)  SRVC 

URITE(6,207)  FLAIdL, IMAT), SLAKIL, IMAT)  SRVC 

....CALCULATE  THE  MEAN  PROJECTION  (OP)  IN  THE  DIRECTION  OF  THE  SOURCE  SRVC 
...(THETA)  OF  ONE  UNIT  LEAF  AREA  WITH  INCLINATION  INCLF .  THE  LEAVES  SRVC 
....AT  THIS  AN6LE  ARE  ASSUMED  TO  BE  AZIMUTHALLY  ISOTROPIC.  SRVC 

DO  330  IANGLE*1 , NSOUR  SRVC 


INCLF«-5.*CEDTR  SRVC 

DO  320  1-1,9  SRVC 

INCLF«INCLF*DELF  SRVC 

320  CALL  COP! INCLF, THETA! I  ANGLE >, OP ( I ) )  SRVC 

C... .CALCULATE  THE  MEAH  PROJECTION  ( OPM )  IN  THE  DIRECTION  OF  THE  SOURCE  SRVC 
C.... (THETA)  OF  ONE  UNIT  LEAF  AREA  AVERA6ED  OVER  THE  CANOPY  LEAF  ANGLE  SRVC 

C... .DENSITY  FUNCTION  F.  SRVC 

CALL  COPN(F,OP,OPM(IANGLE>)  SRVC 

C.... CALCULATE  THE  PROBABILITY  OF  A  HIT  <PHIT>  FOR  A  LIGHT  RAY  UITH  SRVC 

C - SOURCE  DIRECTION  THETA.  SRVC 

CALL  PDENS! IL, IHAT , IAN6LE,QPM( I AN6LE) )  SRVC 

WRITE(4,235)  OP, OPM! I ANGLE) ,PHIT ( IL,IMAT , I ANGLE)  SRVC 

330  CONTINUE  SRVC 

351  CONTINUE 

350  CONTINUE  SRVC 

URITE(6,228)  SRVC 

C.... REFLECTANCE  AND  TRANSMISSION  VECTORS  ARE  READ  FOR  EACH  CANOPY  SRVC 

C - CONSTITUENT.  IN  ADDITION  REFLECTANCE  VECTORS  ARE  READ  FOR  THE  SOIL  SRVC 

C.... BACKGROUND  AND  THE  MEASURED  CANOPY.  THE  MEAN  VECTOR  AND  COVARIANCE  SRVC 
C. . . .AND  CORRELATION  MATRICES  ARE  CALCULATED  AS  UELL  AS  THE  SQUARE-ROOT  SRVC 
C.... MATRIX  WHICH  IS  SUBSEQUENTLY  USED  FOR  MULTIVARIATE  NORMAL  SRVC 

C... .STOCHASTIC  VECTOR  SAMPLING.  SRVC 

C  SRVC 

C.... WAVELENGTHS  TO  BE  SIMULATED  SRVC 

REAB( IFILE, 1 01 )  <XLAM< I ) ,1*1 ,NLAM)  SRVC 

WRITE < & , 20 1 >  (XLAM( I ) ,1*1 ,NLAM)  SRVC 

C.... CONSTITUENT  OPTICAL  VECTORS  SRVC 


C. . .READ  NUMBER  OF  CONSTITUENT  OPTICAL  VECTORS  WHICH  EQUALS  2»MTYPE 
C  *  NUMBER  OF  LAYERS 

READ  < IF  1 LE , 1 05 )  NOP 
105  FORMAT  (110) 

READ ( IFILE, 104) (DATAID(I),I*1,7) 

WRITER, 5)  <DATAID!I),I*1,7) 

READt IFILE, 101)  (CANRH! J),J«1 ,NLAM) 

WRITE<6,2?4) 

WRITE<4,203)  (CANRM! J) , J»1 ,NLAM> 

READ(IFILE,101 )  <SKTIH< J) , J-1 ,NLAM) 

WRITE ( 6 , 295 ) 

WRITE! 4,203)  (SKYIHU) ,  J«1  ,NLAM) 

READ! IFILE, 101)  (DIFIM! J) ,  J*1 ,NIAH) 

URITE<4,2?4) 

WRITE! 6,203 )  (DIFIM! J) , J-1 ,NLAM) 

READ! IFILE, 101 )  (RG( J) , J*1 ,NLAH> 

WRITE! 6,297) 

WRITE(4,203)  (RG( J), J»1 ,NLAM) 
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DO  11  NL*1,NOP 

READ(IFIL£,104)  NULAY ,NTYP, RORT,  ( DAT  A  ID  <  I ) , 1*1 ,7) 
WRITE(6,202)(DATAID<I),I*1 ,7) f HULA Y , HT YP , RORT 
READ<  IFILE, 101 )  (XHU(I,NL),I*1 ,NLAN) 

INL(NULAY,HTYP,RORT)=NL 

WRITE (6, 204)  (XNU( I ,NL) , 1*1 ,NLAH) 

11  CONTINUE 
WRITE  (6,210) 

C .  BIO  LOOP  . 

IST0P*0 

DO  (0  J*1,NLAH 

SOURCE < 1 , J ) »<  SKY IM<J)-DIFIH(J))/<SKYIM(J) ) 

DO  40  1*1 , NBANDS 

40  S0URCE(I+1,J)*(DIFIM(J)*XK(1))/(SKYIH(J>> 

WRITE (6,209) 

DO  43  I*1,NS0UR 

43  WRITE<4,203)  ( SOURCE < 1 , J) , J*t , MLAM) 

DO  7000  ISAHP*I,NSAHP 
DO  4000  ITRIAL*1 , NTRIAL 

C.... COMPUTE  PROPORTION  OF  IRRADIANCE  WHICH  IS  DIRECT  AND  PROPORTION 
C... .WHICH  IS  DIFFUSE. 

C... .POPULATE  FIRST  (TOP)  DOWN  DWELL  LAYER  (DR)  WITH  INCIDENT  DIRECT  AND 
....DIFFUSE  LIGHT.  DOWN  DWELL  RADIATION  FLUX  (DR)  IS  INDEXED  FROM  1  TO 
. . . .NLAY  IN  A  DOWN  GOING  SEQUENCE.  UPWARD  DWELL  RADIATION  FLUX  (UR) 
....15  INDEXED  FROM  1  TO  NLAY+l  IN  UPWARD  GOING  SEQUENCE.  THAT  IS  FOR 
....FOR  UR,  LAYER  1  IS  THE  LAYER  IMMEDIATELY  ABOVE  THE  BACKGROUND.  THE 
....FLUX  IN  LAYER  NLAY+1  IS  THAT  WHICH  ESCAPES  THE  CANOPY  AND  TOGETHER 
....WITH  THE  INCIDENT  FLUX  DETERMINES  THE  CANOPY  REFLECTANCE. 

DO  8  K»1 , 17 
DO  8  J*1,10 
DO  8  1*1,4 

UR(I,J,K)*0. 

DR(I,J,K)*0. 

IGQD< I , J)*0. 

IGOU( I , J)*0. 

8  CONTINUE 
DO  1003  J*1 ,NSOUR 
DO  1003  K*1,NLAN 
1003  DR< 1 , J,K)=SQURCE( J,K) 

C. . . .SET  FLUX  LEVEL  INDICATORS  (DOWNWARD) 

CALL  £THRES(NLAY,NS0UR,-1 ) 

C . FAST  LOOP  TRACES  LIGHT  ATTENUATION  THROUGH  CANOPY . 

C. . . .FLUX  PASSING  THROUGH  LAYERS  IN  A  DOWNWARD  DIRECTION 
2000  CONTINUE 

DO  2600  IL*1,NLAY 


SRVC 

SRVC 

SRVC 

SRVC 


SRVC 

SRVC 

SRVC 

SRVC 

SRVC 

SRVC 

SRVC 

SRVC 

SRVC 

SRVC 

SRVC 


SRVC 

SRVC 

SRVC 

SRVC 

SRVC 

SRVC 

SRVC 

SRVC 

SRVC 
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DO  2500  JS0R=1 , NSOUR  SRVC 

C. . . .CHECK  FLUX  LEVEL  ’NDICATOR  SRVC 

IFdGODdL,  JSOR)  .EQ.O. )  00  TO  2500  SRVC 

C....DID  LIGHT  STRIKE  LEAF  SRVC 

CALL  PGAPdL, JSOR f -1  ,IHIT,MTYPE)  SRVC 

IF<  IHIT  .EQ.O >  GO  TO  2200  SRVC 

DO  2100  IPHIPa1 , 1 8  SRVC 

C.. ..DIRECTION  COSINES  OF  SOURCE  SECTOR  (LVLH)  SRVC 

SXL  *  XS< JSOR, IPHIP)  SRVC 

SYL  a  YStJSOR, IPHIP)  SRVC 

SZL  *  ZS( JSOR)  SRVC 

CALL  LAMBTN< IL, JS0R,NTYPE,-1 ,NSOUR)  SRVC 

2100  CONTINUE  SRVC 

GO  TO  2400  SRVC 

C....GAP  ENCOUNTERED  IN  DOUNUARD  PATH  SRVC 

2200  DO  2250  KL=1 ,NLAH  SRVC 

2250  DR(IL  +  1 , JSOR, KL)aDR(lL+1 , JSOR,KL)+DR<IL, JSOR, XL)  SRVC 

2400  CALL  SETZUL, JS0R,-1)  SRVC 

2500  CONTINUE  SRVC 

CALL  ETHRES(NLAY, NSOUR, -1)  SRVC 

2600  CONTINUE  SRVC 

C.... BACKGROUND  REACHED  -  REFLECTS  LAMBERTIAN  SRVC 

DO  3600  JSOR® 1 , NSOUR  SRVC 

DO  3400  JJa2, NSOUR  SRVC 

IL  -  NLAY  ♦  1  SRVC 

DO  3400  KLa1 , NLAH  SRVC 

UR(  1 , JJ,KL  >aUR<  1 ,  JJ,XL)  +R6(KL)*DRdL,  JSOR,KL)*XK<  JJ-1 ) 

3400  TABS0(4,KL)aTABSQ(4,KL)*< 1 ,-R6(KL) >*DR( IL , JSOR,KL)*XK( JJ-1 ) 

CALL  SETZ < ML AY  + 1 , JSCR,-1 )  SRVC 

3600  CONTINUE  SRVC 

CALL  ETHRES(NLAY, NSOUR, +1)  SRVC 

C. . . .FLUX  PASSING  THROUGH  LAYERS  IN  AN  UPUARD  DIRECTION  SRVC 

DO  4600  IL*1 ,NLAY  SRVC 

DO  4500  JS0R=2, NSOUR  SRVC 

C.... CHECK  FLUX  LEVEL  INDICATOR  SRVC 

IFIIGOUdL, JSOR). EQ.O)  SO  TO  4500  SRVC 

Cr • . .DID  LIGHT  STRIKE  LEAF  SRVC 

CALL  PGAPdL, JSOR,  +  1, IHIT, MTYPE)  SRVC 

IF ( IHIT. EQ.O)  GO  TO  4200  SRVC 

DO  4100  IPHIP-1',  1 8  SRVC 

C.... DIRECTION  COSINES  OF  SOURCE  SECTOR  (LVLH)  SRVC 

SXL  »  XS< JSOR, IPHIP)  SRVC 

SYL  *  YS< JSOR, IPHIP)  SRVC 

SZL  «  ZS(JSOR)  SRVC 

CALL  LAHITNdL, JSOR, MTYPE, +1 , NSOUR)  SRVC 


4100  CONTINUE  SRVC 

GO  TO  4400  SRVC 

C....GAP  ENCOUNTERED  IN  UPWARD  PATH  SRVC 

4200  DO  4250  XL=1,NLAN  SRVC 

4250  URUL  +  1, JSOR, XL)*URULH, JSOR, XL)+UR(IL,JSQR, XL)  SRVC 

4400  CALL  SETZ(IL,JS0R,+1>  SRVC 

4500  CONTINUE  SRVC 

CALL  ETHRES(NLAY,NSQUR, H )  SRVC 

4000  CONTINUE  SRVC 

CALL  ETHRES(NLAY,HS0UR,-1 )  SRVC 

CALL  ETHR£S<NLAY,NS0UR,+1)  SRVC 

C.... RECYCLE  THROUGH  LAYERS  UNTIL  FLUX  EXHAUSTED  SRVC 

DO  5000  IL=1,NLAY  SRVC 

DO  5000  JS0R=2,NS0UR  SRVC 

IF  ( IGOU< IL, JSOR) .NE.O)  GO  TO  2000  SRVC 

5000  CONTINUE  SRVC 

DO  5001  ILa2,NLAYP1  SRVC 

DO  5001  JS0R*1,NS0UR  SRVC 

IFIIGODUL, JSOR). NE.O)  SO  TO  2000  SRVC 

5001  CONTINUE  SRVC 

C _ FLUX  EXHAUSTED  IN  ALL  SOURCES — COMPUTE  REFLECTANCE  FOR  THIS  TRIAL  SRVC 

DO  5200  JSQR=2,NS0UR  SRVC 

DO  5200  KLa1 ,NLAH  SRVC 

RIT ( JSOR, XL) aUR(NLAY+1, JSOR ,KL)/XKUS0R-1 ) 

5200  RITBARC JS0R,XL)=RIT3AR< JSQR,XL)+RIT( JSOR, XL)  SRVC 

URITE<6,283)  ISAHP,  ITRIAL 

DO  5300  JS0Ra2,NS0UR  SRVC 

ZDEG*105-10*JS0R  SRVC 

5300  URITE<6,284)  ZDEG, ( RIT  ( JSOR, XL) ,KLa1 , NLAH)  SRVC 

6000  CONTINUE  SRVC 

C... .TRIALS  COHPLETE  FOR  THIS  SAMPLE  POINT  SRVC 

FTRIAL»NTRIAL  SRVC 

6200  DO  6300  JSQRa2,NSQUR  SRVC 

DO  6300  KL=1 f NLAH  SRVC 

6300  RIT3AR1 JSOR,KL)*RITBAR<  JSOR,XL)/FTRIAL  SRVC 

URITE<6,286)  ISAHP  SRVC 

DO  6400  JS0Ra2,NS0UR  SRVC 

ZDEG»105-10*JS0R  SRVC 

6400  URITE<6,284)  ZDEG, ( RIT BAR ( JSOR, XL ) , KL- 1 f NLAH)  SRVC 

DO  6600  JS0R32,NS0UR  SRVC 

DO  6500  KL»1,NLAH  SRVC 

RBAR(JSOR,KL)«RBARUSOR,KL)+RITBAR(JSOR,KL)  SRVC 

DO  6500  KLL»1,NLAH  SRVC 

6500  COV< JSOR, XL, KLL) *CQV< JSOR, XL, XLL)+RIT8AR< JSOR, XL )*RITBAR< JSOR, KLL)  SRVC 

DO  6600  KL*1 ,NLAH  SRVC 
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4600  RITBARUS0R,KL)*0. 

I F ( I STOP . E Q . 1 )  GO  TO  7100 
7000  CONTINUE 

FSANP*NSAMP 
GO  TO  7150 
7100  FSAHP*1SAMP 

C. . . .ALL  SAMPLE  POINTS  ESTIMATED 
7150  DO  7200  JS0R»2,NS0UR 
DO  7200  KL*1,NLAM 

7200  RBAR( JSOR,KL)*RBAR( JSOR,KL)/FSAMP 
DO  7900  JSQR*2,NS0UR 
ZD£G*105-10*JS0R 
IF (FSAMP.LE.1 . )  GO  TO  7600 
DO  7400  1*1 , NLAM 
DO  7300  J*1 ,NLAM 

7300  C0VUS0R,I,J)a(C0V(JS0RfI,J)-FSAMP»R3AR(JS0R,I)«ftBARUS0RTJ)> 

1/ (F3AMP-1 . ) 

7400  SIG<I)=SORT<COV<JSORfI,I>) 

DO  7500  1*1 , NLAM 
DO  7500  J*1 , NLAM 

7500  COR< I , J)=COM( JSOR, I , J)/ (SIG< I )*SIG( J) ) 

7600  URITE(6,287)  ZDEG, (RBAR( JSOR,KL) ,KL=1 ,NLAM) 

URITE<6,293)  (SIG(KL) ,KL*t rNLAM) 

IF1FSAHP.LE.1.)  GO  TO  7900 
URITE<6,288) 

DO  7700  1*1  , NLAM 

7700  UR1TE(6,289)  ( CQV< JSOR, I , J) , J*1 ,NLAM) 

URITE<6,291 ) 

DO  7800  1*1 , NLAM 

7800  URITE<4,289)  ( COR ( I ,  J) ,  J*1 ,NLAN) 

7900  CONTINUE 

DO  7213  IK*1 ,4 
DO  7213  KL*1,NLAM 

TABSO  < IK,KL)*TABSO( IK,Kl)/(FSAHP*FTRIAL) 

7213  CONTINUE 

DO  7215  1*1,4 
JJ*I 

WRITE  (6,7214)  JJ,  (TABSOU,  J) ,  J*1  ,NLAM) 

7214  FORMAT  (  *  THE  LAYER  IS*, 12,*  THE  ABSORPTIONS  ARE  *, 4(F8.5, 1 X) ) 
7213  CONTINUE 

IF< IFILE.EQ.5)  GO  TO  8000 
STOP 

C . DATA  FORMATS . 

100  F0RHAT<8A10,/,4X,I3,7X,J4,7X,2I2,6X,F6.2,7X,F7.2,5X,F7.2,8X,I2,/, 
15X,I2,7X,I1,7X,IS,9X,IS,8X,I3) 


101 

FORMAT (8F1 O.S) 

SRVC 

102 

FORMAT ( 1 1 0 , 7  A 1 0 ) 

SRVC 

103 

104 

FORMAT  <  8E 10.4) 

F0RHAT<7A10) 

SRVC 

5  FORMAT <#0*f7A10) 

106  FORMAT  4  311 , 7 A 1 0  > 

200  FORMAT<  *1  * , 43X v  *SOLAR  RADI AT  I  ON/ VEGETATION  CANOPY  REFLECTANCE  NODE  SRVC 


1L*,//,44X,*INPUT  DATA*,//, IX, 8A 10,/,  SRVC 

2*  JULIAN  DAY  *,13,*,  YEAR  *,I4,«,  TINE  *,212,*  HOURS*,/,  SRVC 

3*  LATITUDE  »  *,F4.2,*  DEGREES,  LONGITUDE  *  *,F7.2,*  DEGREES*,/,  SRVC 
4*  SOLAR  DECLINATION  =  *,F6.2,*  DEGREES*,/,  SRVC 

S*  BAND  UIDTH  OF  DIFFUSE  VECTORS  -  *,F5. 1  ,*  DEGREES*,/,  SRVC 

6*  NUMBER  OF  UAVELENGTH  BANDS  SIMULATED  *,I2,/,  SRVC 

7*  NUMBER  OF  CANOPY  CONSTITUENTS  *,I1,/,  SRVC 

3*  K  DIGIT  ODD  NO.  TO  INITIALIZE  RANDOM  SEQUENCE  =  *,I5,/,  SRVC 

?*  NSAMP  S*,I3,/,  SRVC 

A*  NTRIAL  =  *,IS,/,  SRVC 

B*  NLAY  *  *,I1,  SRVC 

C)  SRVC 

201  FORMAT(*OUAVELEN6THS  SIMULATED*, /,*0*,F7. 4, 16F8. 4)  SRVC 

202  F0RMAT(*0*,7A10/*  *,*NUMBER  OF  LAYERS  =  *,I1/*  *, 

1 *HATERIAL  TYPE  *  *,I1/*  *,*R  OR  T  *, 1 1 ) 

203  FORMAT (•  *,F7.4,16F3.4)  SRVC 

204  FORMAT < *0  MEAN*,/,3X,10E12.4>  SRVC 

205  FORMAT (*0  COVARIANCE  MATRIX*)  SRVC 

206  FORMAT <*ORANDOM  VECTOR  GENERATED  FROM  THE  *,7A10,/,(*  *,10E12.4))  SRVC 

207  FORMAT(*OLAI  *  *,F4.2,4X,*S  *  *,F4.2)  SRVC 

208  FORHAT<*ODIFFUSE  VECTOR  COEFFICIENTS*,/,  SRVC 

1?<*  K  *) ,/, <?F8.4) )  SRVC 

20?  FORMAT < *OIRRABI ANCE  SOURCE  VECTORS*)  SRVC 

210  FORMAT  < 1H1 )  SRVC 

211  FORMAT (*0  CORRELATION  MATRIX*)  SRVC 

212  F0RHAT(*0DM1  »  *,9F8.4)  SRVC 

221  FORMAT < *OTHRESD  *  *,10F3.4/*  THRESU  *  *,10F8.4)  SRVC 

222  FORMAT < *ODIRECTION  COSINES  OF  SUN  *,3F9.4)  SRVC 

223  FORNAT<*OTHETA  =*, 1 0F8.4)  SRVC 

227  FORMAK///*  *,25< 1H. )  ,2X,*CAN0PY  GEOMETRY*, 2X, 25( 1H. )//)  SRVC 

228  FORMATS*  *,25<IH.))  SRVC 

230  FORMAT < *OLEAF  ANGLE  COMPUTATIONS  -  IL  *  *,I1,  SRVC 

1*  1MAT  *  *,11,*  NANG  «  «,12,/,*  XLFA,TLFA*,  SRVC 

1 /, <2X, 1 6F8.3) )  SRVC 

231  FORMAT < *ONANGLE ( IL,1MAT)  «  *,I2)  SRVC 

232  F0RMAT(*0  FLA  »*,10F3.3)  SRVC 

233  F0RMAT(*0  F  »*,10F8.3)  SRVC 

233  F0RMAT(*0  OP  «*,?F8.3,3X,*0PM  *  *,F8.3,3X,*PHIT  *  *,F8.3)  SRVC 


251  FORMAT (8X, 1 0E1 2.4 )  SRVC 

252  FORMAT <*0REFER  *  *,8E13.4)  SRVC 

283  FORMAT < ‘©REFLECTANCE  FOR  SAMPLE*, 13,*  TRIAL*, 13, 5X,  SRVC 

1 ‘COMPUTATION  TIME  UAS*,F3.1,»  SECONDS.*)  SRVC 

284  FORMAT <*  2  »*,I3,*  DEG*,3X,tQF7.3>  SRVC 

285  FORMAT(*OCAUTION... .SAMPLE*, 13,*  CONTAINS  ONLY*, 13,*  TRIALS.*)  SRVC 

286  FORMAT <*0*,75(tH. )/*  MEAN  REFLECTANCE  FOR  SAMPLE*, 13)  SRVC 

287  F0RMAT(*06RAND  MEAN  FOR  Z  »*,I3,*  DEGREES. *,3X,10F7. 3)  SRVC 

288  FQRHAT(*0  COVARIANCE  MATRIX*)  SRVC 

289  F0RMAT(7X,t0F12.8)  SRVC 

291  F0RMAT(»0  CORRELATION  MATRIX*)  SRVC 

292  FORMAT  < 1 X , 1 20 ( 1 H— ) )  SRVC 


293  FORNAT(*OSTDEV*, 1 X, 10F12.8) 

294  FORMAT(*OMEASURED  REFLECTANCE*) 

295  FORMAT (*06L08AL  IRRADIANCE*) 

296  FORHAT(*ODIFFUSE  IRRADIANCE*) 

297  FORMAT(*OSOIL  REFLECTANCE*) 


END  SRVC 

*DECX  LAMBIN 

SUBROUTINE  LAMBTN< IL, JSOR,MTYPE, IDIR,NSOUR)  LAHBTN 

C . FOR  A  GIVEN  FLUX  SOURCE  THIS  PROGRAM  CALLS  THE  APPROPRIATE  LAHBTN 

C . PR06RAHS  TO  DETERMINE  LEAF  ORIENTATION  AND  OPTICAL  PROPERTIES  LANBTN 

C . AND  UPDATES  THE  DIFFUSE  SOURCES  UITH  SCATTERED  FLUX.  LAHBTN 

C  SXL,  SYL,  SZL  LAHBTN 

C  JSOR  LAHBTN 

C  LXS,  LYS,  LZS  LAHBTN 

C  IDIR  LAHBTN 

C  NLAM  LAHBTN 

C  DR( I, J,K>  LAHBTN 

C  UR< I , J,K)  LAHBTN 

C  MTYPE  LAHBTN 

C  IL  LAHBTN 

C  NSOUR  LAHBTN 

C  R,T 

C  ZENITH  LAHBTN 

C  OUTPUT  VARIABLES  LAHBTN 

C  DR(I,J,K)  LAHBTN 

C  UR< I , J,K)  LAHBTN 

C  LAHBTN 

COMMON  BUH01 17) ,XK<9) ,SXL,SYL,SZL,XLF, YLF,ZLF  LAHBTN 

COMMON/C 1 /DUN  1 <7) ,NLAM,DUM2 (26) ,CE1PI,DUN3(24) ,LXS,LYS,LZS  LAHBTN 

C0MHQN/C2/DUM4<31 ),R(17),T<17),RG<17) ,DUM5< 1 97)  ,ZENITH( 10 ) 

C0MM0N/C6/DR <  4 , 1 0 , 1 7 ) , UR  <  4 , 1 0 , 1 7 )  LAHBTN 

C0MH0N/CNAT/MTP<3) ,NLAY,0PH( 10)  LAHBTN 

C 
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COMMON  / AB3/7 ABSO (4,17) 

DIMENSION  H(17),P7RP<2,17) 

REAL  LXS,LYS,LZS  LAMBTN 

DATA  PI02/1 .570796327/  LAMBTN 

C. . . .SET  DIRECTION  COSINES  OF  SOURCE  LAMBTN 

XL*SXL  LAMBTN 

YLaSYL  LAMBTN 

ZL=SZL  LAMBTN 

IF(JSOR.NE.t)  SO  TO  100  LAMBTN 

XL=LXS  LAMBTN 

YL*LYS  LAMBTN 

ZLaLZS  LAMBTN 

C.... RANDOM  LEAF  ORIENTATION,  DIRECTION  COSINES  OF  NORMAL,  AND  LAM3TN 

C. . . .LEAF  OPTICAL  PROPERTIES  LAMBTN 

100  IFUDIR.EQ.-1)  IXL=IL  LAMBTN 

IF ( IDIR.E3. I )  IXL*NLAY-IL+1  LAMBTN 

CALL  LAN6L£(IXL,MTYPE,THETAL,PHIL)  LAMBTN 

C _ SET  SIDE  OF  LEAF  UHICH  LIGHT  STRIKES.  I SIDE= 1  < TOP ) ,  -1  (BOTTOM).  LAMBTN 

ISIDE=-IDIR  LAMBTN 

DOT»XL*XLF+YL*YLF+ZL»ZLF  LAMBTN 

IFIDOT.LT.O.)  ISIDE*IDIR  LAMBTN 

COSLS*ABS(DQT)  LAMBTN 

IF< IDIR.EQ. 1 )  GO  TO  5  LAMBTN 

DO  A  XL “ 1 , NLAM  LAMBTN 

4  H<KL)»  DR(IL,JS0R,KL)/18. 

GO  TO  9  LAMBTN 

3  DO  7  KL*1,NLAM  LAMBTN 

7  H(KL)*  UR(IL, JS0R,KL)/18. 

9  CONTINUE  LAMBTN 

C...SET  OPTICAL  PROPERTIES  FOR  MTYPE, LAYER, REFLECT.  AND  TRANS. 

C.... UPDATE  DIFFUSE  SOURCES  WITH  SCATTERED  RADIATION  FLUX  LAMBTN 

DO  50  JJS0R*2,NS0UR  LAMBTN 

IF(ISIDE.EQ.-I)  CALL  BFLUX< THETAL, ZENITH ( JJSOR) ,H,T,R,NLAH,PTRP)  LAMBTN 
IF( IS1DE.EQ. 1 )  CALL  BFLUX < THETAL , ZEM I T H ( JJSOR) ,H,R,T, NLAM, PTRP)  LAMBTN 
DO  50  KL»1,NLAH  LAMBTN 

IF ( IDIR.EQ. 1 )  GO  TO  45  LAMBTN 

DR< IL  +  1,  JJSOR, KL)»DR(IL*1 f JJSOR, KL)+PTRP(2,KL)  LAMBTN 

UR(NLAY+2'IL, JJSOR, KL)=UR( NLA Y+2" I L, JJSOR, XL) +PTRP( 1 ,KL)  LAMBTN 

GO  TO  SO  LAMBTN 

45  DR<NLAY+2-IL, JJSOR, XL) =DR(NLAY+2-IL, JJSOR, XL) +PTRP<2, XL)  LAMBTN 

UR( IL+1, JJSOR, XL )*URUL+1, JJSOR, KL )+PTRP(1, XL)  LAMBTN 

50  CONTINUE  LAMBTN 

DO  S3  KL*1,NLAM 

TABS04 IL,KL)*H(KL)  *< 1 .-(R(KL) +T(XL) ) )+TABSO< IL,KL) 

33  CONTINUE 
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RETURN 

END 

'DECK  BFIUX 

SUBROUTINE  8FLUX(TAfTRP,HfR,T,NLAM,PTRP> 

GIVEN  THE  IRRABIANCE  H  OF  A  LEAF  INCLINED  AT  TA  THIS  PROGRAM 
DETERMINES  THE  FLUX  REFLECTED  AND  TRANSMITTED  INTO  A  SOURCE 
BAND  UHOSE  ZENITH  ANGLE  IS  TRP. 

DIMENSION  PTRP<2,17),H<17),R<!7>,T(17) 

DATA  PI/3.1 A15*2054/,P102/1 .570794327/ 
F1<X,r>»C0S<TA)'<SIN<X)'*2-SIN<r>"2> 
F2(X)»AC0S(-1/(TAN<TA)*TAN(X) ) ) 

F3(X,T,Z)'2.*SIN<TA)*SIN<X>*<DEL+.25*<SIN(2.*Y)-SIN( 2.*Z) ) )/PI 
DEL*. 087264443 
T1-TRP-DEL 
T2*TRP+D£L 

IF1TA.LE.PIQ2-T2)  GO  TO  10 
IF1TA.GE.PIQ2-T1)  GO  TO  20 
GO  TO  30 
C....CASE  1 

10  XF1»FMT2,T1) 

DO  15  KL»1,NLAH 
PTRP< 1 ,KL)*R(KL)'H<XL)'XF1 
PTRP(2,KL)*T(XL)*H(KL)'XF1 
15  CONTINUE 
RETURN 
C. . . .CASE  2 

20  XF1»F1 <T2,Tl ) 

IF1TA.LE. 1.5533)  GO  TO  21 
PRP*PI02 

60  TO  22 

21  PRP-F21TRP) 

22  XF3«F3(PRP,T1,T2) 

DO  25  KL*1,NLAM 

PTRP< 1  ,KL)«H(KL)*(R(KL)+T (KL) )'XF3* 
1(R(KL)*H<KL)'PRP-T(KU»H<KU*PI+T(KL)*H(KL)*PRP)*XF1/PI 
PTRP(2,KL)«H(KL)*(T(KL)*R<KL) )*XF3* 

1(T<KL)*H(KL)'PRP-R(KL)*H(XL>'PI*-R(XL)*H(KL)*PRP)*XF1/P1 
25  CONTINUE 
RETURN 
C....CASE  3 
30  TB*PI02-TA 
XFl'Fl <  TB  r  T 1 ) 

DO  35  KL»1 ,NLAM 
PTRP( 1 ,KL)»R(KL)*H<KL)*XF1 
35  PTRP(2,KL)»T <KL)*H(KL)*XF1 
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IF< TB+T2.LE.3. 1 06 )  60  TO  36  BFLUX 

PRP-PI02  BFLUX 

60  TO  37  BFLUX 

36  PRP*F2( (TB+T2)/2. )  BFLUX 

37  XF1«F1(T2,TB>  BFLUX 

DEL«< <TRP*TA)/2.)-. 74176493  BFLUX 

XF3*F3(PRP,TB,T2)  BFLUX 

DO  40  KL«1,NLAM  BFLUX 

PTRP( 1  ,KL)=PTRP<  1 ,KL)+H(KL)*(R(KL)+T (KL) )*XF3+  BFLUX 

1 (R<KL)*H(KL)*PRP-T (KL)*H(KL)*PI+T(KL>*H<KL)*PRP)*XF1 /PI  BFLUX 

PTRP(2,KL)SPTRP(2,KL)+H(KL)*(T (KL)+R(KL>)*XF3*  BFLUX 

1 <T(XL)*H(KL)*PRP-R<XL)*H<KL)*PI+R(KL)*H(KL)*PRP)*XF1/PI  BFLUX 

40  CONTINUE  BFLUX 

RETURN  BFLUX 

ENB  BFLUX 

♦DECK  LANGLE 

SUBROUTINE  LANGLE* IL,HTYPE, THETAL, PHIL)  LANGLE 

——THIS  PR06RAN  SELECTS  A  RANDON  LEAF  INCLINATION  (THETAL)  AND  AZIMUTH  LANGLE 
(PHIL)  AND  THEN  COHPUTES  ITS  DIRECTION  COSINES  XLF ,  YLF ,  AND  2LF.  LAN6LE 

THE  INTERMEDIATE  PARAMETERS  SINL,  COSL,  SIN P,  AND  COSP  ARE  ALSO  LANGLE 

OUTPUT.  RANDON  LEAF  REFLECTANCE  AND  TRANSMITTANCE  VECTORS  ARE  ALSO  LAN6LE 
SELECTED.  LANGLE 

LANGLE 

INPUT  LAN6LE 

II  LANGLE 

NTYPE  LANGLE 

NANGLE  LANGLE 

OUTPUT  LANGLE 

THETAL  LANGLE 

PHIL  LANGLE 

XLF,  YLF,  ILF  LANGLE 

SINL,  COSL,  SINP,  COSP  LANGLE 

M 

LANGLE 

COHHON/C1 / DUH2  <31 ) ,CERTD,DUH?<3) ,CE2P1  LANGLE 

C0HH0N/C4/NAM6LE(3,3) ,FLA(3,3, 1 0) , SLA I (3, 3) ,FLAI(3,3),PHIT(3,3,10)  LANGLE 
CQMM0N/C8/SINL, COSL, SINP, COSP  LANGLE 

COMMON  DUM3<29) ,XLF, YLF,2LF  LANGLE 

COMMON  /KIM/  I NL (3,3,2) 

C - DETERMINE  RANDOM  LEAF  ORIENTATION.  LANGLE 

FM*NAN6LE ( IL ,HT YPE )  LANGLE 

XT»RANF(0.)  LANGLE 

XI«1  .♦(FM-1 . )*XT  LANGLE 

IX*XI  LANGLE 

IF ( IX. EQ. NANGLE (IL,MT YPE ) )  IX-IX-T  LANGLE 


n  o  o 


IXP1*1X+1 


THETAL=FLA< IL,MTYPE, IX )+.3*<FLA< IL,MTYPE, IXP1 »-FLA< IL ,HTYPE , IX > > 
PH1L=CE2PI*RA#F<0.) 

C - THE7AL,  PHIL  ARE  LEAF  INCLINATION  AND  AZIMUTH,  RESPECTIVELY. 

1  CONTINUE 

SINL=SIN(THETAL) 

COSL»COS( THETAL) 

SINP*SIN(PHIL) 

COSP=COS(PHIL) 

C  — -COMPUTE  LEAF  NORMAL  DIRECTION  COSINES 
XLF»-SINL«COSP 
YLF«-SINL*SINP 
ZLF*COSL 

C - SELECT  RANDOM  LEAF  REFLECTANCE  AND  TRANSMITTANCE  VECTORS. 

CALL  OPTICAUMTYPE, IL) 

RETURN 

END 

♦  DECK  COP 

SUBROUTINE  COP< ALPHA, BETA, OP) 

C 

C....THIS  PROGRAM  CALCULATES  THE  MEAN  PROJECTION  OF  A  UNIT  LEAF  AREA  IN 

C _ THE  DIRECTION  OF  THE  SOURCE.  THE  LEAF  IS  INCLINED  AT  AN  ANGLE 

C... .ALPHA  AND  IS  ASSUMED  TO  BE  AZIMUTHALLY  ISOTROPIC.  THE  SOURCE 
C... .DIRECTION  IS  AT  AN  AZIMUTH  OF  ZERO  AND  AN  INCLINATION  OF  BETA. 

C 


COMMON/Ct/DUNl (33) ,CEP102 
OPaCOS( ALPHA ) ♦SIN <  BETA) 

IF( ALPHA. LE. BETA)  RETURN 

C....THETAO  IS  THE  LEAF  AZIMUTH  ANGIE  AT  UHICH  OP  3EC0MES  NEGATIVE  AND 
C....IS  IN  THE  FIRST  QUADRANT.  THE  FUNCTION  OP  IS  SYMMETRIC  AND  HENCE 
C....IS  AVERA6ED  OVER  LEAF  AZIMUTH  AN6LES  OF  0  TO  PI  RADIANS. 
THETAO»ACOS( TAM (BETA) /TAN (ALPHA) ) 

TANTO»TAN(TH£TAO) 

OP»OP*(1.*<TAMTO-THETAO)/CEPIQ2) 

RETURN 

END 

♦DECK  COPM 

SUBROUTINE  COPM(G, OP, OPM) 

C. . . .THIS  PROGRAM  CALCULATES  THE  MEAN  PROJECTION  OF  A  UNIT  LEAF  AREA  IN 

C _ THE  DIRECTION  OF  THE  SOURCE  (OPM)  FOR  THE  SIMULATED  CANOPY.  THE 

C.... LEAVES  OF  THE  CANOPY  ARE  ASSUMED  TO  BE  AZIMUTHALLY  ISOTROPIC.  THE 
....OP  FUNCTION  USED  IN  THE  CALCULATION  HAS  BEEN  PREVIOUSLY  DETERMINED 
....FOR  A  GIVEN  SOURCE  DIRECTION  FOR  LEAF  INCLINATION  ANGLES  OF 
....3,  15,  ...,  85  DEGREES.  G  IS  THE  LEAF  INCLINATION  ANGLE  DENSITY 
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COPN 

DIMENSION  OP  (?) f6<9) 

COPN 

OPNsO. 

COPN 

DO  1  1*1,9 

COPN 

1 

OPH»OPM+OP(I)*G(I) 

COPN 

RETURN 

COPN 

END 

COPN 

♦  DECK 

PDENS 

SUBROUTINE  PDENS dL, NT YPE, I ANGLE, QPN) 

PDENS 

C— — 

THIS  PROGRAM  COMPUTES  THE  PROBABILITY  THAT  LIGHT 

AT 

INCIDENT  ANGLE 

PDENS 

c 

THETA(IANGLE)  INTERACTS  UITH  MATERIAL  TYPE  NTYPE 

UITHIN 

CANOPY 

PDENS 

c 

LAYER  IL. 

PDENS 

c 

PDENS 

c 

INPUT 

PDENS 

c 

IL 

PDENS 

c 

NTYPE 

PDENS 

c 

IANGLE 

PDENS 

c 

OPH 

PDENS 

c 

SLAI 

PDENS 

c 

FLAI 

PDENS 

c 

THETA 

PDENS 

c 

OUTPUT 

PDENS 

c 

PHIT 

PDENS 

c 

PDENS 

C0NN0N/C2/DUN(28?) , THETA( 10) 

C0NH0N/C4/NANGLE(3,3) , FLA (3, I, 1 0) , SLAI (3,3) , FLAI (3 

,3),! 

PHIT (3,3,10) 

PDENS 

ARG*1 .-(SLAI < IL, NTYPE )*OPN/S IN ( THETA ( I ANGi £) ) ) 

PDENS 

IF  (ARG.LE.O.)  GO  TO  1 

PDENS 

PO*ARG^*(FLAI ( IL, NTYPE) /SLAI 1 IL, NTYPE) ) 

PDENS 

GO  TO  2 

PDENS 

1 

PO  *  0. 

PDENS 

URITE(6,100)  IANGLE 

PDENS 

too 

FORMAT  ( 1  HO,  ♦  PO  SET  TO  ZERQM5) 

PDENS 

2 

CONTINUE 

PDENS 

PHITdL, NTYPE, 1ANGLEM.-P0 

PDENS 

RETURN 

PDENS 

END 

PDENS 

♦DECK 

PGAP 

SUBROUTINE  PGAPdL, IANGLE, ID1R,IHIT, NTYPE) 

PGAP 

c— - 

THIS  PROGRAM  DETERMINES  IF  AN  INTERACTION  IS  BEING 

MADE 

IN  LAYER  IL 

PGAP 

C 

AND  SETS  THE  MATERIAL  TYPE  OF  LAYER  IL. 

PGAP 

C 

PGAP 

C 

INPUT 

PGAP 

C 

IL 

PSA  P 

c 

IANGLE 

PGAP 

noon 


C  1 0 1 R  PGAP 

C  NLAY  PGAP 

C  NTP  P GAP 

C  PM I T  P6AP 

C  OUTPUT  POAP 

C  IHIT  PGAP 

C  NTYPE  PGAP 

C  PGAP 

C0HN0N/C4/NANGLE<3,3),FLA<3,3, 10)  ,SLAI (3,3) ,FLAI (3,3) ,PHIT< 3,3,10)  PGAP 
COMMON/ CN AT / NTP(3) ,NL AY  P GAP 

REAL  PHITN 

IF(IDIR.LT.O)  GO  TO  10  PGAP 

ILAYER»NLAY+1 -IL  PGAP 

GO  TO  20  PGAP 

10  ILAYER°IL  PGAP 

C _ MTP ( ILAYER )  GIVES  THE  LAST  NTYPE  WITHIN  A  LAYER  WHICH  CONTAINS  THE  COMBINED 

C  NTYPE  DISTRIBUTION. 

20  NTYPE*MTP(ILAYER)  PGAP 

IHIT=0  PGAP 

TEST=RANF(0. )  PGAP 

IF(PHIT< I LAYER, NTYPE, I ANGLE ) .LT.TEST )  GO  TO  30  PGAP 

IHIT=1  PGAP 

IF  (NTYPE. EQ.  1 )  GO  TO  30 

C...A  HIT  HAS  BEEN  RECORDED  -  NOU  WHAT  DID  IT  HIT. 

C.. .NORMALIZE  THE  2  MATERIAL  DISTIBUTION  TO  1.0 


PHITN  «  PHIT ( IL, 1 , I ANGLE) /(PHI T( IL,  1 ,IANGLE)+PHIT(lL,2,lANGL£)) 
TEST  *RANF<0. ) 

IF  (PHITN  .LT.TEST)  GO  TO  40 
NTYPE  »  1 
GO  TO  30 
40  MTYPE  ■  2 


30  RETURN  P GAP 

END  PGAP 

•DECK  ETHRES 

SUBROUTINE  ETHRES (NL A Y, NSOUR , IDIR)  ETHRES 

—THIS  PROGRAM  DETERMINES  (FOR  EACH  LAYER  AND  FOR  ALL  LIGHT  SOURCE  ETHRES 

DIRECTIONS)  IF  THE  SOURCE  FLUX  IS  ABOVE  THRESHOLD  REQUIREMENTS  IN  ETHRES 

THE  DIRECTION  INDICATED  BY  IDIR.  INDICATORS  IGOD  OR  IGOU  ARE  SET  ETHRES 
ACCORDINGLY.  ETHRES 

C  ETHRES 

C  INPUT  ETHRES 

C  NLAY  ETHRES 

C  NSOUR  ETHRES 

C  IDIR  ETHRES 

C  NLAM  ETHRES 


ooooooo  #  —  o  o  n  o  o  o  o  n 


:  DR  ETHRES 

UR  ETHRES 

THRES  ETHRES 

OUTPUT  ETHRES 

IGOD  ETHRES 

I60U  ETHRES 

ETHRES 

COHNON/C1 /DUN0<7) ,NLAM  ETHRES 

C0HH0H/CA/!)R(4,10f17>,UR<4,10f17),THRESD<10),IG0D(4,t0),IG0U(4,10)  ETHRES 
1 ,  THRESUUO)  ETHRES 

—-DOWNWARD  FLUX  ETHRES 

IF ( ID1R.6T .0)  GO  TO  10  ETHRES 

NLAYER3NLAY+1  ETHRES 

DO  2  I»t ,HLAYER  ETHRES 

DO  2  J=1 ,  NSOUR  ETHRES 

IGOD< I , J)=0  ETHRES 

DO  1  K*1 f  NLAH  ETHRES 

IF(DR<I,J,K).LT.THRESD< J))  GO  TO  1  ETHRES 

IGOD(I  ,J)»1  ETHRES 

GO  TO  2  ETHRES 

1  CONTINUE  ETHRES 

2  CONTINUE  ETHRES 

RETURN  ETHRES 

- UPWARD  FLUX  ETHRES 

0  CONTINUE  ETHRES 

DO  4  I31,NLAY  ETHRES 

DO  4  J*2f NSOUR  ETHRES 

IG0U(I,J)«0  ETHRES 

DO  3  K*1 ,  NLAH  ETHRES 

IF(UR(I,J,K).LT.THRESU<J))  GO  TO  3  ETHRES 

IG0U(I,J)*1  ETHRES 

GO  TO  4  ETHRES 

3  CONTINUE  ETHRES 

4  CONTINUE  ETHRES 

RETURN  ETHRES 

END  ETHRES 

DECK  SETZ 

SUBROUTINE  S£TZ( IL, IANGLE, IDIR)  SETZ 

- THIS  PR06RAN  SETS  THE  FLUX  <AND  ITS  APPROPRIATE  INDICATORS)  IN  THE  SETZ 

IDIR  DIRECTION  AT  ANGLE  THETA(IANGLE)  IN  LAYER  IL  TO  ZERO.  SETZ 

SETZ 

INPUT  SETZ 

IL  SETZ 

IANGLE  SETZ 

IDIR  SETZ 


/ 
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C  NLAM  SETZ 

OUTPUT  SETZ 

DR  SETZ 

UR  SETZ 

IfiOD  SETZ 

I60U  SETZ 

SETZ 

C0MN0N/C1/DUN1<7),NLAN  SETZ 

CQHH0N/Ci/DR<4, 10, 1 7) ,UR<4, 1 0, 1 7) ,THRES< 1 0) , IG0B< 4, 1 0) , IG0U( 4 , 10)  SETZ 
IFdDIR.EQ.1 )  GO  TO  10  SETZ 

— DOWNWARD  FLUX  SETZ 

DO  1  K>1,NLAM  SETZ 

1  DRdL,IAN6LE,K)*0.  SETZ 

1 GOD ( ILv1ANGL£>30  SETZ 

RETURN  SETZ 

— — UPWARD  FLUX  SETZ 

0  CONTINUE  SETZ 

DO  2  K»1,NLAN  SETZ 

2  URdL,  IAN6L£,K)S0.  SETZ 

IGOUdL,IAN6LE)=0  SETZ 

RETURN  SETZ 

END  SETZ 

DECK  OPTICAL 

SUBROUTINE  OPTICAL  (MTYPE, IL> 

——THIS  PROGRAH  SELECTS  RANDON  LEAF  REFLECTANCE  AND  TRANSMITTANCE  OPTICAL 

VECTORS  FOR  MATERIAL  TYPE  MTYPE.  OPTICAL 

OPTICAL 

INPUT  OPTICAL 

MTYPE  OPTICAL 

NVEC  OPTICAL 

C  OPTICAL 

XMU  OPTICAL 

OUTPUT  OPTICAL 

R*T 

OPTICAL 

C0MMQN/L1/DATAID(7,?),XMU<17,?>,C(17,17,?),NVEC<?)  OPTICAL 


C0MM0M/C2/CANRH(17),SKYIM<17),DIFIN(t7>rR(17),T(17>, 

1 RG< 17) , XL AM  < 1 7 ) ,SOURCE( 1 0, 1 7) r  THE  T  A ( 1 0 ) 

COMMON  /KIM/  INL<3,3,2) 

C.. .SELECT  APPROPRIATE  OPTICAL  VECTOR  GIVEN  MTYPE, IL, AND  R  OR  T  VECTOR. 

I»  INL(IL, MTYPE, 1) 

J«I  +  1  OPTICAL 

It  CALL  UTIL(XMU( 1 , I ) ,R) 

13  CALL  UTIL(XMU(1,J),T) 

RETURN  OPTICAL 
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OPTICAL 


END 

♦DECK  NRH 
♦DECK  MATSQR 
♦DECK  81 DATA 

SLOCK  DATA  BDAT 

COMM0N/C1/DUH(3O) ,CEDTR,CERTD,CEMTR,C£PIQ2,CE1PI ,CE2PI  3DAT 

DATA  C£DTR,CERTD,CEMTR/. 017453293, 57. 2957795, .00029088821/  3DAT 

DATA  C£PI02,CE1PI,C£2PI/1. 37079632, 3. 14159263, 6. 23318530/  BDAT 

END  3DAT 

♦DECK  TBLR 

SUBROUTINE  TBLR(H,  X,  Y,  XX,  Z)  TBLR 

C  T3LR 

C. . . .THIS  PROGRAM  FINDS  THE  INTEGRAL  Z(X)  OF  THE  FUNCTION  Y(X)  FROH  X<1)  TBLR 

C _ TO  X12N-1)  USING  SIMPSONS  RULE  THE  INTEGRAL  Z(X)  IS  NORMALIZED  TO  TBLR 

C _ 1.0  AT  X ( 2M- 1 ) .  THE  TABLE  OF  Z  VERSUS  X  IS  THEN  INVERTED  TO  DETER-  T3LR 

C....NINE  X  AS  A  FUNCTION  OF  Z  AT  M  REGULARLY  SPACED  POINTS  ALONG  Z.  T3LR 

C  TBLR 

C  INPUT  VARIABLES  TBLR 

C  M  =  DESIRED  NUMBER  OF  REGULARLY  SPACED  POINTS  ALONG  Z  TBLR 

C  X  =  SPECIFIED  AT  2M-1  POINTS  TBLR 

C  Y  »  SPECIFIED  AT  2M-1  POINTS  TBLR 

C  OUTPUT  VARIABLES  TBLR 

C  XX  »  THE  TABLE  OF  X  VALUES  FOR  M  REGULARLY  SPACED  POINTS  T3LR 

C  (H-1  INTERVALS)  ALONG  Z.  TBLR 

C  Z  a  THE  NORMALIZED  INTEGRAL  OF  Y  AT  X(1),  X<3),  ...,  X(2M-1).  TBLR 

C  TBLR 

DIMENSION  xm>,Y<1?>,Z()0),XI<10),XX<10>  TBLR 

C... .SIMPSONS  RULE  INTEGRATION  TBLR 

10  Z<1>  »  0.0  TBLR 

DX  «  X< 2)  -  X(1)  TBLR 

20  DO  50  J  =»  2,M  TBLR 

JO  =»  2^J  -  3  TBLR 

30  J1  *  2^J  -  2  TBLR 

J2  *  2*J  -  1  TBLR 

AO  ZU>  *  Z(J  -  1)  ♦  DX*<Y<  JO)  M.*Y<  J1 )  ♦  Y  <  J2  >  )/3.0  TBLR 

50  XI ( J)  *  X < J 2 )  TBLR 

XI < 1 >*X< 1 )  TBLR 

C... .NORMALIZE  INTEGRAL  Z(X)  TBLR 

60  DO  70  J  «  1,M  TBLR 

70  Z(J)  »  Z(J)/Z<M)  TBLR 

C. . . .FIND  X  AT  M  REGULARLY  SPACED  POINTS  ALONG  Z.  TBLR 

XX(1)  »  III)  TBLR 

EM  »  N  -  1  TBLR 

F  «  1.0/EN  TBLR 

JS*2  TBLR 
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)0  DO  120  K  >  2,N 

2T  -  K  -  t 
2T  -  ZT*F 

70  DO  110  J  »JS,H 

I F < 2 < J >  -  ZT)  110,  100,  100 
100  G  *  (ZT  -  Z(J  -  1>>  /(Z( J)  -  Z ( J  -  1)> 

XX ( K )  =  XI < J  -  1)  ♦  0*<XI(J)  -  XI ( J  -  1)> 

GO  TO  115 
110  CONTINUE 
115  JS*J 
120  CONTINUE 
RETURN 
END 

SUSROUTINE  SUN 

——THIS  PROGRAM  CALCULATES  THE  POSITION  OF  THE  SUN 

INPUT 

TINE 

6LAT 

DEC 

OUTPUT 

SINLAT,  COSLAT 
SINDEC,  COSDEC 
COSH 

SINZ,  COSZ 
SINAZ,  COSAZ 
LXS,  LYS,  LZS 

TINE  OF  SIHULATION  (HOURS) 

GLAT  IS  SITE  GEOGRAPHICAL  LATITUDE 

6L0N6  IS  SITE  LONGITUDE 

DEC  IS  SOLAR  DECLINATION 

H  IS  SOLAR  HOUR  ANGLE 

COSZ  IS  COSINE  OF  SOLAR  ZENITH  AN6LE 

COSAZ  IS  COSINE  OF  SOLAR  AZIMUTH 

LXS,  LYS,  LZS  ARE  SOLAR  DIRECTION  COSINES 


C0HN0N/C1/DAY, YEAR, TINE, GLAT, GLONG, DEC, DUN( 24), 
1CEDTR,CERTD,CEHTR,DUM2(I?>, 

2SINLAT, COSLAT, SINDEC, COSDEC, COSH, SINZ, COSZ, SINAZ, COSAZ, LXS, LYS, LZS 
REAL  LXS, LYS, LZS 
H*ABS<  <<12.-TIME)*15.)*CEDTR> 

SINLAT*SIN(GLAT) 


/ 


c 


COSLAT-COS<  GLAT )  SUN 

SIND£C*SIN<D£C)  SUN 

COSDECaCOS( DEC)  SUN 

COSH»COS<H>  SUN 

C0S2«SINLAT*S1NBEC*C0SLAT*C0SDEC*C0SH  SUN 

SIN2«SQRT(1.-C0SZ*C0SZ)  SUN 

C0SA2«(SIND£C-SINLAT*C0SZ)/<CQSLAT»SIHZ>  SUN 

SINA2*SQRT ( I .-C0SAZ*C0SA2)  SUN 

LXS*S1N2*C0SAZ  SUN 

LYS«SIN2*SINAZ  SUN 

LZS»COSZ  SUN 

RETURN  SUN 

END  SUN 

C 

SUBROUTINE  UTIL<A,B)  UTIL 

C.....SET  VECTOR  B  *  VECTOR  A  UTIL 

C0HN0N/C1/DUN(7) ,NLAM  UTIL 

OIHENSION  A < 1 7 ) ,3(17)  UTIL 

DO  1  I»1,NLAH  UTIL 

1  B(I)*A(I)  UTIL 

RETURN  UTIL 

END  UTIL 

C 

SUBROUTINE  FUN(A,B)  FUN 

RETURN  FUN 

END  FUN 


A53 
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APPENDIX  B 


GEOMETRICAL  MATRICES  FOR  THEORETICAL  CANOPIES 


B1 


Planoohile  Canopy  Geometry 


FOLIAGE  ANGLE  DISTRIBUTION 


INCLINATION  ANGLE 

9R0IAIIUTV  OF  OCCURRENCE 
LAYER  1  LAYER  2  LAYER  3 

0.0 

.105 

.105 

.105 

5.0 

.104 

.104 

.104 

10.0 

.102 

.102 

.102 

15.0 

.098 

.098 

.098 

20.0 

.093 

.093 

.093 

25.0 

.086 

.086 

.086 

10.0 

.079 

.079 

.079 

35.0 

.071 

.071 

.071 

40.0 

.062 

.062 

.062 

45.0 

.053 

.053 

.053 

50.0 

.043 

.043 

.043 

55.0 

.035 

.035 

.035 

60.0 

.026 

.026 

.026 

65.0 

.019 

.019 

.019 

70.0 

.012 

.012 

.012 

75.0 

.007 

.007 

.007 

80.0 

.003 

.003 

.003 

85.0 

.001 

.001 

.001 

90.0 

0.000 

0.000 

0.000 

B2 


CANOPY  GEOMETRY  INPUT  DATA  FOR  PLAN0PHILE,LAI-1 


LEAF  AREA  INDEX 
LAYER  1  LAYER  2  LAYER  3 
.25  .50  .25 


CANOPY  DENSITY  PARAMETERS 
LAYER  I  LAYER  2  LAYER  3 
.10  .10  .10 

GEOMETRICAL  VIEW  ANGLE  FACTOR  MATRIX 
ZENITH  AN6LE 


5.0 

15.0 

25.0 

35.0 

45.0 

55.0 

65.0 

75.0 

85.0 

LAYER  1 

.5006 

.2652 

.2094 

.1913 

.1846 

.1820 

.1812 

.1810 

.1810 

LAYER  2 

.3546 

.3380 

.2964 

.2798 

.2732 

.2707 

.2699 

.2697 

.2696 

LAYER  3 

.0476 

.1052 

.1035 

.1012 

.1001 

.0996 

.0995 

.0994 

.0994 

GROUND 

.0373 

.2916 

.3907 

.4277 

.4422 

.4477 

.4495 

.4499 

.4499 

LONG  HAVE  TRANSFER  MATRIX 


FROM 


TO 

SKY 

LAYER  1 

LAYER  2 

LAYER  3 

GROUND 

LAYER  1 

.8695 

.2540 

.3385 

.1087 

.4222 

LAYER  2 

.5973 

.1663 

.4657 

.3663 

.5973 

LAYER  3 

.4222 

.1087 

.3385 

.2540 

.8695 

B3 
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CANOPY  8C0HETRV  INPUT  DATA  FOR  PLANOPHILE,LAI»4 

LEAF  AREA  INDEX 
LAYER  t  LAYER  2  LAYER  3 
1.00  2.00  1.00 

CANOPY  DENSITY  PARAMETERS 
LAYER  1  LAYER  2  LAYER  3 
.10  .10  .10 


GEOMETRICAL  VIEW  ANGLE  FACTOR  MATRIX 
ZENITH  ANGLE 


3.0 

13.0 

23.0  33.0 

43.0 

55.0  63.0  75.0 

85.0 

LAYER  1 

.9627 

.7084 

.6093  .3723 

.3578 

.5523  .5505  .5301 

.5501 

LAYER  2 

.0372 

.2668 

.3311  .3493 

.3537 

.3379  .3587  .3588 

.3589 

LAYER  3  0 

.0000 

.0176 

.0363  .0448 

.0482 

.0495  .0500  .0501 

.0501 

GROUND  0 

.0000 

.0072 

.0233  .0335 

.0382 

.0402  .0408  .0410 

.0410 

LONG  UAVE 

TRANSFER 

MATRIX 

FROM 

TO 

SKY 

LAYER  1  1 

LAYER  2 

LAYER  3  GROUND 

LAYER 

1 

.5973 

.7983 

.4944 

.0588  .0441 

LAYER 

2 

.1578 

.2196 

1  .2381 

.2196  .1578 

LAYER 

3 

.0441 

.0588 

.4944 

.7983  .5973 

B4 
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CANOPY  GEOMETRY  INPUT  DATA  FOR  PLANOPHILE, LAI-7 


LEAF  AREA  INDEX 
LAYER  t  LAYER  2  LAYER  3 
1 .73  3.30  1.73 


CANOPY  DENSITY  PARAMETERS 
LAYER  1  LAYER  2  LAYER  3 
.10  .10  .10 


GEOMETRICAL 

VIEU  ANGLE 

FACTOR  MATRIX 

ZENITH  AN8LE 

5.0 

15.0 

23.0 

35.0  45.0 

35.0 

65.0  75.0 

LAYER  1 

.9968 

.8843 

.8069  . 

7738  .7602 

.7550 

.7333  .7529 

LAYER  2 

.0032 

.1141 

.1859  . 

2146  .2260 

.2303 

.2317  .2321 

LAYER  3  0 

.0000 

.0014 

.0058  . 

0090  .0105 

.0111 

.0113  .0114 

GROUND  0 

.0000 

.0002 

.0014  . 

0026  .0033 

.0036 

.0037  .0037 

LONG  WAVE  TRANSFER  MATRIX 

FROM 

TO 

SKY 

LAYER 

1  LAYER  2 

LAYER 

3  GROUND 

LAYER 

1 

.4222 

1.1484 

.4011 

.0162 

.0049 

LAYER 

2 

.0441 

.1516 

1.6014 

.1516 

.0441 

LAYER 

3 

.0049 

.0162 

.4011 

1.1484 

.4222 

B5 
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Erectophile  Canopy  Geometry 

FOLIAGE  ANGLE  DISTRIBUTION 


PROBABILITY  OF  OCCURRENCE 


INCLINATION  ANGLE 

LAYER  1 

LAYER  2 

LAYER  3 

0.0 

0.000 

0.000 

0.000 

5.0 

.001 

.001 

.001 

10.0 

.003 

.003 

.003 

15.0 

.007 

.007 

.007 

20.0 

.012 

.012 

.012 

25. 0 

.019 

.019 

.019 

30.0 

.024 

.024 

.024 

35.0 

.035 

.035 

.035 

40.0 

.043 

.043 

.043 

45.0 

.053 

.053 

.053 

50.0 

.042 

.042 

.042 

55. 0 

.071 

.071 

.071 

40.0 

.079 

.879 

.079 

45.0 

.084 

.084 

.084 

70.0 

.093 

.093 

.093 

75.0 

.098 

.098 

.098 

80.0 

.102 

.102 

.102 

85.0 

.104 

.104 

.104 

90.0 

.105 

.105 

.105 

CANOPY  GEOMETRY  INPUT  DATA  FOR  ERECTOPHILE,LAIB1 


LEAF  AREA  INDEX 
LAYER  1  LAYER  2  LAYER  3 
.25  .50  .25 


CANOPY  DENSITY  PARAMETERS 
LAYER  1  LAYER  2  LAYER  3 
.10  .10  .10 


6E0METRICAL  VIEW  ANGLE  FACTOR  MATRIX 


5.0 

15.0 

25.0 

ZENITH 

35.0 

AN6LE 

45.0 

55.0 

65.0 

75.0 

85.0 

LAYER  1 

.9056 

.4333 

.2770 

.2036 

.1620 

.1365 

.1204 

.1109 

.1064 

LAYER  2 

.0936 

.3847 

.3451 

.2912 

.2495 

.2196 

.1991 

.1863 

.1801 

LAYER  3 

.0000 

.0769 

.1047 

.1026 

.0953 

.0879 

.0020 

.0780 

.0759 

GROUND 

.0001 

.1031 

.2732 

.4024 

.4931 

.5561 

.5985 

.6248 

.6375 

LONG  NAVE  TRANSFER  MATRIX 


FROM 


TO 

SKY 

LAYER  1 

LATER  2 

LAYER  3 

GROUND 

LAYER  1 

.7591 

.4716 

.3378 

.0864 

.3350 

LAYER  2 

.4764 

.1540 

.7292 

.1540 

.4764 

LAYER  3 

.3350 

.0864 

.3378 

.4716 

.7591 

B7 


CANOPY  GEOMETRY  INPUT  DATA  FOR  ERECTOPHILE, LAI-4 


LEAF  AREA  INDEX 
LAYER  f  LAYER  2  LAYER  3 
1.00  2.00  1.00 

CANOPY  DENSITY  PARAMETERS 
LAYER  1  LAYER  2  LAYER  3 
.10  .10  .10 


GEOMETRICAL  VIEW  ANGLE  FACTOR  MATRIX 


3.0 

13.0 

23.0 

ZENITH 

33.0 

ANGLE 

45.0 

53.0 

65.0 

75.0 

85.0 

LAYER  1 

.9999 

.0969 

.7268 

.5976 

.5069 

.4439 

.4015 

.3752 

.3623 

LAYER  2 

.0001 

.1020 

.2320 

.3372 

.3732 

.3841 

.3841 

.3809 

.3784 

LAYER  3 

0.0000 

.0010 

.0148 

.0389 

.0608 

.0763 

.0861 

.0915 

.0939 

GROUND 

0.0000 

.0001 

.0036 

.0262 

.0591 

.0936 

.1283 

.1524 

.1652 

LONG  HAVE  TRANSFER  MATRIX 


FROM 


TO 

SKY 

LAYER  1 

LAYER  2 

LAYER  3 

GROUND 

LAYER  1 

.4764 

1.0372 

.3794 

.0467 

.0502 

LAYER  2 

.1377 

.1627 

1.3891 

.1627 

.1377 

LAYER  3 

.0502 

.0467 

.3794 

1.0372 

.4764 

B8 
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CANOPY  GEOMETRY  INPUT  IATA  FOR  IRfCTQPHILf ,LAI»7 


LEAF  AREA  INDEX 
LAYER  1  LAYER  2  LAYER  3 
1.73  3.30  1.75 

CANOPY  DENSITY  PARAMETERS 
LAYER  1  LAYER  2  LAYER  3 
.10  .10  .10 


6E0METR1CAL  VIEIi 

ANSI! 

FACTSR  MATRIX 

ZENITH 

ANBLE 

3.0 

15.0 

25.0  35.0 

45.0 

55.0  05.0 

75.0 

83.0 

LAYER  1  1.0000 

.9812 

.8907  .7967 

.7098 

.0419  .5928 

.5009 

.5451 

LAYER  2  0.0000 

.0188 

.1022  .1949 

.2057 

.3122  .3397 

.3544 

.3007 

LAYER  3  0.0000 

0.0000 

.0010  .0007 

.0173 

.0295  .0400 

.0475 

.05f3 

8R0UND  0.0000 

0.0000 

.0001  .0017 

.0071 

.0104  .0275 

.0372 

.0428 

LONG  WAVE 

TRANSFER  MATRIX 

FROM 

TO 

SKY 

LAYER  1  LAYER  2 

LAYER  3  6R0UND 

LAYER  1 

.3350 

1.3200 

.3057 

.0188 

.0105 

LAYER  2 

.0502 

.  1  f  5 1  1 

.0592 

.1151 

.0502 

LAYER  3 

.0105 

.0188 

.3057 

1.3200 

.3350 

B9 
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Plagiophi le  Canopy  Geometry 


FOLIAGE 


INCLINATION  AN6LE 

0.0 

5.0 

10.0 

15.0 

20.0 

25.0 

30.0 

35.0 

40.0 

45.0 

50.0 

55.0 

60.0 

65.0 

70.0 

75.0 

80.0 

85.0 

90.0 


ANGLE  DISTRIBUTION 

PROBABILITY  OF  OCCURRENCE 


LAYER  1 

LAYER  2 

LAYER  3 

0.000 

0.000 

0.000 

.003 

.003 

.003 

.013 

.013 

.013 

.028 

.028 

.028 

.046 

.046 

.046 

.065 

.065 

.065 

.083 

.083 

.083 

.098 

.098 

.098 

.108 

.108 

.108 

.111 

.111 

.111 

.108 

.108 

.108 

.098 

.098 

.098 

.083 

.083 

.083 

.065 

.065 

.065 

.046 

.046 

.046 

.028 

.028 

.028 

.013 

.013 

.013 

.003 

.003 

.003 

0.000 

0.000 

0.000 

CANOPY  GEOMETRY  INPUT  DATA  FOR  PLA6I0PHILE, LAI-1 


LEAF  AREA  INDEX 
LAYER  1  LAYER  2  LAYER  3 
.23  .30  .25 


CANOPY  DENSITY  PARAMETERS 
LAYER  t  LAYER  2  LAYER  3 
.10  .10  .10 


GEOMETRICAL  VIEW  ANGLE  FACTOR  MATRIX 


3.0 

13.0 

25.0 

ZENITH 

33.0 

ANGLE 

45.0 

55.0 

65.0 

75.0 

85.0 

LAYER  1 

.S005 

.3622 

.2461 

.2003 

.1001 

.1716 

.1684 

.1676 

.1676 

LAYER  2 

.1916 

.37S4 

.3254 

•  2SG3 

.2660 

.2599 

.2563 

.2557 

.2556 

LAYER  3 

.0064 

.0940 

.1033 

.1024 

.0993 

.0975 

.0969 

.0967 

.0967 

GROUND 

.0016 

.1654 

.3230 

.4090 

.4318 

.4710 

.4782 

.4800 

.4802 

LQN6  WAVE 

TRANSFER 

MATRIX 

FROM 

TO 

SKY 

LAYER  1 

LAYER  2 

LAYER  3 

6R0UND 

LAYER  1 

.8124 

.3680 

.3525 

.0976 

.3622 

LAYER  2 

.5212 

.1665 

.6173 

.1665 

.5212 

LAYER  3 

.3622 

.0976 

.3525 

.3680 

.8124 

CANOPY  GEOMETRY  INPUT  DATA  FOR  PLA6I0PHILE, LAI-4 


LEAF  AREA  INDEX 


LAYER  1 

LAYER  2 

LAYER  3 

1.00 

2.00 

1.00 

CANOPY 

DENSITY  PARAMETERS 

LAYER  1 

LAYER  2 

LAYER  3 

.10 

.10 

.10 

GEOMETRICAL  VIEU 

AN6LE  FACTOR  MATRIX 

ZENITH 

AN8LE 

5.0 

15.0 

25.0  35.0 

45.0 

55.0  65.0  75.0 

85.0 

LAYER  I  .9984 

.8346  . 

6770  .5910 

.5482 

.5290  .5218  .5200 

.5198 

LAYER  2  .0016 

.160?  . 

2893  .3406 

.3596 

.3665  .3688  .3694 

.3695 

LAYER  3  0.0000 

.0038  . 

0228  .0404 

.0506 

.0553  .0571  .0575 

.0575 

SROUND  0.0000 

.0007  . 

0109  .0280 

.0417 

.0492  .0523  .0531 

.0532 

LONS  HAVE  TRANSFER 

MATRIX 

FROM 

TO 

SKY 

LAYER  1  LAYER  2 

LAYER  3  GROUND 

LAYER  1 

.5212 

.9503 

4337 

.0491  .0384 

LAYER  2 

.1335 

.1892  1. 

3475 

.1892  .1335 

LAYER  3 

.0384 

.0491 

4337 

.9503  .5212 

CANOPY  GEOMETRY  INPUT  DATA  FOR  PLAGIOPHILE,LAI<*7 


LEAF  AREA  INDEX 
LAYER  1  LAYER  2  LAYER  3 
1.73  3.50  1.75 


CANOPY  DENSITY  PARAMETERS 
LAYER  1  LAYER  2  LAYER  3 
.10  .10  .10 


6E0METRICAL  VIEU 

AN6LE 

FACTOR  MATRIX 

ZENITH 

ANGLE 

5.0 

15.0 

25.0  35.0 

45.0 

55.0  65.0 

75.0 

85.0 

LAYER  1  1.0000 

.9571 

.8616  .7908 

.7510 

.7322  .7250 

.7232 

.7230 

LAYER  2  0.0000 

.0428 

.1357  .2000 

.2335 

.2486  .2542 

.2556 

.2557 

LAYER  3  0.0000 

.0001 

.0023  .0072 

.0116 

.0141  .0151 

.0153 

.0154 

6R0UND  0.0000 

0.0000 

.0004  .0019 

.0038 

.0051  .0057 

.0059 

.0059 

LONS  HAVE 

TRANSFER  MATRIX 

FROM 

TO 

SKY 

LAYER  t  LAYER  2 

LAYER  3  GROUND 

LAYER  1 

.3622 

1.2683 

.3433 

.0142 

0047 

LAYER  2 

.0384 

.1271  1 

.6617 

.1271 

0384 

LAYER  3 

.0047 

.0142 

.3433 

1.2683 

3622 

B13 
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Extremophile  Canopy  Geometry 


FOLIAGE  ANGLE  DISTRIBUTION 


PROBABILITY  OF  OCCURRENCE 


INCLINATION  AN6LE 

LAYER  1 

LAYER  2 

LAYER  3 

0.0 

.100 

.100 

.100 

S.O 

.097 

.097 

.097 

10.0 

.088 

.088 

.088 

ts.o 

.075 

.075 

.075 

20.0 

.059 

.059 

.059 

23.0 

.041 

.041 

.041 

30.0 

.025 

.025 

.025 

33.0 

.012 

.012 

.012 

40.0 

.003 

.003 

.003 

45.0 

0.000 

0.000 

0.000 

50.0 

.003 

.003 

.003 

55.0 

.012 

.012 

.012 

60.0 

.025 

.025 

.025 

65.0 

.041 

.041 

.041 

70.0 

.059 

.059 

.059 

75.0 

.075 

.075 

.075 

80.0 

.088 

.088 

.088 

85.0 

.097 

.097 

.097 

00.0 

.100 

.100 

.100 

B14 
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CANOPY  GEOMETRY  INPUT  DATA  FOR  EXTREM0PHILE,LAI»1 


LEAF  AREA  INDEX 
LAYER  1  LAYER  2  LAYER  3 
.23  .30  .23 


CANOPY  DENSITY  PARAMETERS 
LAYER  t  LAYER  2  LAYER  3 
.10  .10  .10 


GEOMETRICAL  VIEU  AN6LE 

FACTOR  MATRIX 

ZENITH  ANGLE 

3.0 

15.0 

25.0  33.0  45.0 

55.0 

65.0  75.0 

85.0 

LAYER  1 

.7427 

.3437 

.2414  .1949  .1672 

.1484 

.1354  .1270 

.1228 

LAYER  2 

.2403 

.3736 

.3221  .2832  .2552 

.2340 

.2183  .2077 

.2023 

LAYER  3 

.0126 

.0972 

.1054  .1017  .0966 

.0917 

.0875  .0845 

.0829 

GROUND 

.0044 

.1035 

.3311  .4202  .4811 

.5259 

.5588  .5808 

.5920 

LONG  NAVE  TRANSFER  MATRIX 

FROM 

TO 

SKY 

LAYER  1  LAYER  2 

LAYER 

3  GROUND 

LAYER 

1 

.8350 

.3202  .3383 

.0971 

.3995 

LAYER 

2 

.5572 

.1618  .5521 

.1618 

.5572 

LAYER 

3 

.3995 

.097)  .3383 

.3202 

.8350 

CANOPY  GEOMETRY  INPUT  DATA  FOR  EXTREMOPHILE, LAI-4 


LEAF  AREA  INDEX 


LAYER  1 

LAYER  2 

LAYER  3 

t.00 

2.00 

1.00 

CANOPY 

DENSITY  PARAMETERS 

LAYER  t 

LAYER  2 

LAYER  3 

.10 

.to 

.10 

GEOMETRICAL  VIEU 

ANGLE 

FACTOR  MATRIX 

ZENITH 

ANGLE 

3.0 

13.0 

25.0  35.0 

45.0 

55.0  65.0 

75.0 

85.0 

LAYER  1 

.9956 

.8145 

.6689  .5798 

.5189 

.4741  .4412 

.4192 

.4080 

LAYER  2 

.0044 

.1792 

.2948  .3460 

.3697 

.3804  .3843 

.3849 

.3845 

LAYER  3  0 

.0000 

.0052 

.0243  .0430 

.0578 

.0690  .0770 

.0821 

.0846 

GROUND  0 

.0000 

.0012 

.0120  .0312 

.0536 

.0765  .0975 

.1138 

.1228 

LONG  HAVE 

TRANSFER  MATRIX 

FROM 

TO 

SKY 

LAYER  1  LAYER  2 

LAYER  3  GROUND 

LAYER 

1 

.5572 

.8756 

.4436 

.0576 

0561 

LAYER 

2 

.1643 

.1956  1 

.2704 

.1956 

,1643 

LAYER 

3 

.0561 

.0576 

.4436 

.8756 

,5572 

B16 


CANOPY  GEOMETRY  INPUT  DATA  FOR  EXTREN0PHILE,LAI*7 


LEAF  AREA  INDEX 
LAYER  1  LAYER  2  LAYER  3 
1.75  3.50  1.75 


CANOPY  DENSITY  PARAMETERS 
LAYER  1  LAYER  2  LAYER  3 
.10  .10  .10 


GEOMETRICAL  VIEM 

AN6LE 

FACTOR  MATRIX 

ZENITH 

ANGLE 

5.0 

15.0 

25.0  35.0 

45.0 

55.0  65.0 

75.0 

85.0 

LAYER  1 

.9999 

.9482 

.856*  .7821 

.7236 

.6768  .6404 

.6152 

.6020 

LAYER  2 

.0001 

.0517 

.1404  .2074 

.2550 

.2891  .3127 

.3274 

.3344 

LAYER  3  0 

.0000 

.0001 

.0020  .0061 

.0111 

.0162  .0207 

.0241 

.0259 

GROUND  0 

.0000 

0.0000 

.0010  .0044 

.0103 

.0179  .0262 

.0334 

.0376 

LONG  WAVE 

TRANSFER  MATRIX 

FROM 

TO 

SKY 

LAYER  1  LAYER  2 

LAYER  3  6R0UND 

LAYER 

1 

.3983 

1.1936 

.3674 

.0149 

,0160 

LAYER 

2 

.0558 

.1427  1 

.5931 

.1036 

,0949 

LAYER 

3 

.0125 

.0273 

.5174 

.8756 

,3572 

B17 


t 


% 


Spherical  Canopy  Geometry 

FOLIAGE  ANGLE  DISTRIBUTION 


PROBABILITY  OF  OCCURRENCE 


INCLINATION  ANGLE 

LAYER  1 

LAYER  2 

LAYER  3 

0.0 

.056 

.056 

.056 

5.0 

.056 

.056 

.056 

10.0 

.056 

.056 

.056 

15.0 

.056 

.056 

.056 

20.0 

.056 

.056 

.056 

25.0 

.056 

.056 

.056 

30.0 

.056 

.056 

.056 

35.0 

.056 

.056 

.056 

40.0 

.056 

.056 

.056 

45.0 

.056 

.056 

.056 

50.0 

.056 

.056 

.056 

55.0 

.056 

.056 

.056 

60.0 

.056 

.056 

.056 

65.0 

.056 

.056 

.056 

70.0 

.056 

.056 

.056 

75.0 

.056 

.056 

.056 

80.0 

.056 

.056 

.056 

85.0 

.056 

.056 

.056 

90.0 

.056 

.056 

.056 

CANOPY  GEOMETRY  INPUT  BATA  FOR  SPHERICAL, LAI* t 


LEAF  AREA  INDEX 
LAYER  1  LAYER  2  LAYER  3 
.25  .50  .25 


CANOPY  DENSITY  PARAMETERS 

LAYER  f  LAYER  2  LAYER  3 

.10  .10 

.10 

GEOMETRICAL  VIEU  ANGLE 

FACTOR  MATRIX 

ZENITH  ANGLE 

5.0 

15.0 

25.0  35.0  45.0 

55.0  65.0 

75.0 

LAYER  1 

.7712 

.3525 

.2437  .1975  .1733 

.1594  .1511 

.1464 

LAYER  2 

.2169 

.3760 

.3237  .2856  .2617 

.2467  .2372 

.2316 

LAYER  3 

.0092 

.0957 

.1054  .1021  .0979 

.0947  .0924 

.0911 

6R0UND 

.0027 

.1758 

.3273  .4148  .4670 

.4992  .5192 

.5309 

LONG  HAVE  TRANSFER  MATRIX 

FROM 

TO 

SKY 

LAYER  l  LAYER  2 

LAYER  3  GROUND 

LAYER 

1 

.8244 

.3427  .3455 

.0974 

3814 

LAYER 

2 

.5398 

.1643  .5832 

.1643 

,5398 

LAYER 

3 

.3814 

.0974  .3455 

.3427 

,8244 

B19 


/ 


85.0 

.  1 442 

.2290 

.0904 

.5364 


CANOPY  GEOMETRY  INPUT  DATA  FOR  SPHERICAL ,LA1«4 


LEAF  AREA  INDEX 
LAYER  1  LAYER  2  LAYER  3 
1.00  2.00  1.00 


CANOPY  DENSITY  PARAMETERS 
LAYER  1  LAYER  2  LAYER  3 
.10  .10  .10 

GEOMETRICAL  VIEW  ANGLE  FACTOR  MATRIX 
ZENITH  ANGLE 


5.0 

15.0 

25.0 

35.0 

45.0 

55.0 

45.0 

75.0 

85.0 

LAYER  1 

.9973 

.0242 

.4727 

.5052 

.5330 

.5008 

.4808 

.4491 

.4434 

LAYER  2 

.0027 

.1703 

.2922 

.3435 

.3452 

.3748 

.3792 

.3813 

.3821 

LAYER  3 

0.0000 

.0045 

.0234 

.0418 

.0543 

.0423 

.0473 

.0702 

.0715 

GROUND 

0.0000 

.0010 

.0115 

.0294 

.0474 

.0421 

.0727 

.0794 

.0828 

LONG  HAVE 

!  TRANSFER 

MATRIX 

FROM 

TO 

SKY 

LAYER  1 

LAYER  2 

LAYER  3 

GROUND 

LAYER  1 

.5398 

.9118 

.4395 

.0534 

.0448 

LAYER  2 

.1489 

.1928 

1.3079 

.1928 

.1489 

LAYER  3 

.0448 

.0534 

.4395 

.9118 

.5398 

B20 
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CANOPY  GEOMETRY  INPUT  DATA  FOR  8PHERICAL, LAI-7 


LEAF  AREA  INDEX 
LAYER  1  LAYER  2  LAYER  3 
1.75  3.50  1.75 


CANOPY  DENSITY  PARAMETERS 
LAYER  1  LAYER  2  LAYER  3 
.10  .10  .10 

GEOMETRICAL  VIEW  ANGLE  FACTOR  MATRIX 
ZENITH  ANGLE 


5.0 

15.0 

25.0 

35.0 

45.0 

55.0 

65.0 

75.0 

85.8 

LAYER  1 

1.0000 

.9523 

.8584 

.7856 

.7362 

.7035 

.6824 

.6698 

.6638 

LAYER  2 

0.0000 

.0476 

.1388 

.2046 

.2455 

.2704 

.2855 

.2942 

.2982 

LAYER  3 

0.0000 

.0001 

.0024 

.0077 

.0135 

.0183 

.0219 

.0241 

.0252 

GROUND 

0.0000 

0.0000 

.0004 

.0021 

.0048 

.0077 

.0102 

.0119 

.0128 

B.00  12.00  16.00  20.00  2M.00  28.00  32.00  36.00  HO. 00 

AIR  TEMP  (OEGREES  C) 


GROUND  TEMP  (DEGREES  C) 


Tl.OO  0.20  U. <10  0.60  0.80  J .  00  1.20  1.00  J .  60  1 .  BO 

nELflTIVE  HUMIDITY 


MO. 00  80.00  120.00  160.00  200.00  2M0.00  280.00  320.00  360.00 

W1NU  SPEED  (CM/SEC) 


00*62  00^22  00 '02  CO '  9  T  00 '9 1  OO'fit  CC'2'f5 

(3ayyoiiN23  S33y93Q) 3uruuu3dW3i 


.00  MO. 00  00.00  120.00  160.00  200.00  240.00  280.00  320.00  360.00 

SHORT  WAVE  ABSORPTION  LAYER  I 


tl.00  4 .  Of)  0.00  12.00  16.00  20.00  2U.00  28.00  32.00  36.00 

SHOW  WAVE  ABSORPTION  LATER  2 


00'Ti2- 


00  ‘ZZ 


*T).  88  0.90  0.92  0. 9'1  0.96  0.98  1.00  1.02  1.04 

EMISSIVITY  LAYER  2 


.90  0.92  0. 9*1  0.96  0.98  1.00 

EMISSIVITY  LRYER  3 


SENSITIVITY  PLOT  DF  DRY 


(DE 


-  I  <>.  00  -8.00  -1|.00  0.00 

GROUND  TEMP  (DEGREES  CJ 


40.00  80.00  120.00  160.00  200.00  240.00  280.00  320.00  360.00 

WIND  SPEED  ICH/SEC) 


RE  LOT 


STUMATRL  RFS I  STANCE  CM  I N/CM) 


□o< 


O— Q — B — Q — Q 


LATER 


OF  GROUND 


COEFFICIENT  LAYER 


Oak-Hickory  Daytime  Sensitivity  Plots 


C33 


/ 


.00  >1.00  8.00  12.00  16.00  20.00  24.00  28.00  32.00  36.00 

GROUND  TEMP  (DEGREES  Cl 


ST0MRTRL  RES  I 5TRNCE  CM  I N/CM) 


80.00  120.00  160.00  200.00  280.00  280.00  320.00  360.00 

SHOUT  WAVE  ABSORPTION  LATER  1 


'll. 00  4.00  0.00  12.00  16.00  20.00  24.00  28.00  32.00  36.00 

SHORT  WAVE  ABSORPTION  LATER  2 


U3JUJ1  NOIidyOSflb 


~i).  88  0.90  0.92  0.99  0.96  0.98  1.00  1.02  I.OM  1.06 

EM15S1VITY  OF  GROUND 


ABSORPTION  COEFFICIENT  LATER 


UNCLASSIFIED 

3or3 


COLORADO  STATE  UNXV  FORT  COLLINS  DEPT  OF  FOREST  AND  —ETC  F/G  17/5 
THERMAL  VEGETATION  CANOPY  MOOEL  STUDIES. (U) 

AUG  81  J  A  SMITH*  K  J  R ANSON*  D  NGUYEN  DACW39-77-C-0073 

MES/TR/cL-81-6 


END 


I  ENT 


Oak-Hickory  Nighttime  Sensitivity  Plots 


C50 


•  i  * 


/ 


SENSITIVITY  PLOT  OH  NIGHT 


O 

o 


GObbO  1 1N33  S33y030)  3afUbU3dW3l 


C53 


GROUND  TEMP  (DEGREES  C) 


>10.00  80.00  120.00  160.00  200.00  240.00  280.00  320.00  360.00 

WIND  SPEED  (CM/SEC) 


ST0MATAL  RESISTANCE  (MIN/CM) 


EM1SSIVITT  LATER 


0.80  0.64  0.88  0.92  0.98  1.00  1.04  1.08 

flBSORPT  1  ON  COEFFICIENT  LAYER  2 


APPENDIX  D:  SUPPORTING  VALIDATION  DATA 


D1 


/ 

t 


Cedar  River,  Douqlas-Fir 


CANOPY  GEOMETRY  INPUT  OAT A  FOR  DOUGLAS-FIR 


LEAF 
LAYER  1 
1.50 


AREA  INDEX 
LAYER  2  LAYER  3 
5.30  1.00 


CANOPY  DENSITY  PARAMETERS 
LAYER  1  LAYER  2  LAYER  3 
.10  .10  .10 


FOLIAGE  ANGLE  DISTRIBUTION 


PROBABILITY  OF  OCCURRENCE 


INCLINATION  ANGLE 

LAYER  1 

LAYER  2 

LAYER 

0.0 

.088 

.093 

.089 

5.0 

.078 

.079 

.078 

10.0 

.079 

.080 

.079 

15.0 

.077 

.078 

.078 

20.0 

.084 

.084 

.084 

25.0 

.077 

.077 

.077 

30.0 

.081 

.080 

.080 

35.0 

.059 

.059 

.059 

40.0 

.088 

.087 

.088 

45.0 

.063 

.062 

.063 

50.0 

.062 

.061 

.062 

55.0 

.045 

.043 

.044 

60.0 

.044 

.042 

.043 

65.0 

.029 

.029 

.029 

70.0 

.024 

.024 

.024 

75.0 

.013 

.013 

.013 

80.0 

.007 

.007 

.007 

85.0 

.003 

.003 

.003 

90.0 

0.000 

0.000 

0.000 

D2 


THERHAL  MODEL 

INPUT  DATA  FOR  DOUGLAS-FIR 

GEOMETRICAL  VIEW  ANGLE 

FACTOR  MATRIX 

ZENITH  ANGLE 

5.0 

15.0 

25.0  35.0 

65.0 

55.0  65.0  75.0 

35. < 

LAYER  1 

.9976 

.8703 

.7697  .7216 

.6996 

.6900  .6865  .6856 

.6855 

LAYER  2 

.0026 

.1296 

.2289  .2753 

.2960 

.3067  .3079  .3088 

.3089 

LAYER  3 

.0000 

.000) 

.0009  .0018 

.0025 

.0028  .0029  .0029 

.0029 

GROUND 

.0000 

.0000 

.0006  .0015 

.0022 

.0025  .0027  .0027 

.0027 

LONG  UAVE 

TRANSFER  MATRIX 

FROM 

TO 

SKY 

LAYER  1 

LAYER  2 

LAYER  3  GROUND 

LAYER 

1 

.6599 

1.0731 

.6560 

.0031  .0027 

LAYER 

2 

.0267 

.0661  1 

.8071 

•0507  .0621 

LAYER 

3 

.0022 

.0052 

.5887 

.8008  .5960 

AVERAGE  SHORTUAVE  ABSORPTION  COEFFICIENTS 
LAYER  1  LAYER  2  LAYER  3 
.3890  .0190  ‘.0280 


STOMATAL  RESISTANCE 
.66  <HIN/CM) 


03 


/ 


ENVIRONMENTAL  INPUT  DATA 


£ 


CEDAR  RIVER, 

WASHINGTON 

4  AUGUST  1979 

TIRE 

AIR  TEMP 

GRND  TEMP  UIND  SPEED 

REL  HUH 

6L0BAL  SUR 

(HOURS) 

(DE6  C) 

(DEG  C) 

(H/SEC) 

(U/H**2) 

IOO 

9.8 

12.4 

1.0 

.99 

0.00 

200 

9.3 

12.0 

1.0 

1.00 

0.00 

300 

8.7 

11.5 

.9 

1.00 

0.00 

400 

8.2 

11.1 

1.0 

1.00 

0.00 

500 

7.5 

10.4 

.9 

1.00 

0.00 

400 

7.0 

10.2 

.9 

1.00 

i.ao 

700 

7.2 

9.7 

.7 

1.00 

49.20 

800 

9.2 

9.7 

.4 

.99 

227.30 

900 

11.8 

10.3 

.5 

.94 

445.70 

1000 

13.8 

11.3 

1.5 

.89 

421 .40 

1100 

15.7 

13.0 

1.8 

.87 

772.40 

1200 

14.8 

14.3 

1.9 

.85 

814.50 

1300 

18.3 

15.5 

1.9 

.84 

847.70 

1400 

19.5 

14.2 

2.1 

.83 

484.80 

1500 

20.5 

17.0 

2.2 

.81 

835.70 

1 600 

21.2 

17.9 

2.7 

.81 

770.20 

1700 

21.3 

18.2 

3.0 

.81 

418.40 

1800 

21.3 

18.8 

3.1 

.80 

493.90 

1900 

21.2 

17.7 

2.5 

.78 

289.90 

2000 

29.5 

17.4 

2.0 

.78 

115.40 

2100 

17.4 

14.7 

.8 

.81 

7.50 

2200 

14.4 

15.7 

t.t 

.87 

0.00 

2300 

14.4 

14.8 

1.5 

.88 

0.00 

2400 

14.4 

14.5 

1.5 

.88 

0.00 

D4 


ENVIRONMENTAL  INPUT  DATA 


CEDAR  RIVER,  UASHIHG70N  5  AUGUST  19 79 


TINE  AIR  TENP  GRND  TEMP  U1ND  SPEED  REL  HUN  GLOBAL  SUR 


(HOURS) 

1  AA 

(DEG  C) 

IT  L 

(DEG  C> 

1  J  A 

(N/SEC) 

4  4 

A  A 

(U/N**2) 

A  A  A 

1  V0 

13.6 

1  •  1 

•  72 

0.00 

200 

13.0 

14.2 

.4 

.93 

0.00 

300 

12.7 

14.2 

.5 

.94 

0.00 

400 

12.6 

14.1 

1.1 

.93 

0.00 

500 

12.0 

14.0 

1.2 

.93 

0.00 

600 

11.3 

13.7 

.9 

.94 

1.60 

700 

11.0 

13.4 

.6 

.95 

30.40 

800 

11.2 

13.3 

.6 

.94 

64.50 

900 

11.6 

13.4 

1.2 

.93 

111.90 

1000 

12.2 

13.5 

1.2 

.91 

154.60 

1100 

13.0 

13.7 

1.2 

.90 

228.20 

1200 

14.6 

14.3 

1.3 

.88 

659.60 

1300 

16.7 

15.2 

1.6 

.85 

719.30 

1400 

17.0 

15.8 

2.0 

.84 

370.00 

1500 

17.3 

16.0 

1.3 

.84 

366.60 

1600 

18.6 

16.4 

1.3 

.84 

659.40 

1700 

18.6 

16.8 

1.4 

.83 

399.20 

1800 

19.4 

16.8 

1.8 

.82 

388.00 

1900 

19.4 

16.9 

2.3 

.82 

301.30 

2000 

18.6 

16.6 

1.2 

.83 

110.90 

2100 

15.7 

16.0 

.7 

.87 

9.20 

2200 

13.6 

15.0 

1.0 

.93 

0.00 

2300 

12.5 

14.2 

.9 

.95 

0.00 

2400 

12.4 

13.7 

1.0 

.96 

0.00 

D5 


LAYER  I  LAYER  2  LAYER  3 

3.40  .80  .40 

CANOPY  DENSITY  PARAMETERS 
LAYER  1  LAYER  2  LAYER  3 

.10  .10  .10 


r 

i 


FOLIAGE  ANGLE  DISTRIBUTION 


INCLINATION  ANGLE 

PROBABILITY  OF  OCCURRENCE 
LAYER  1  LAYER  2  LAYER  3 

0.0 

.066 

.117 

.014 

3.0 

.067 

.155 

.233 

10.0 

.084 

.129 

.120 

13.0 

.086 

.177 

.157 

20.0 

.050 

.064 

.053 

25.0 

.098 

.135 

.154 

30.0 

.084 

.081 

.100 

35.0 

.076 

.037 

.047 

40.0 

.063 

.040 

0.000 

45.0 

.087 

.019 

.010 

50.0 

.040 

.015 

0.000 

53.0 

.043 

.019 

0.000 

60.0 

.031 

.007 

0.000 

65.0 

.033 

.002 

0.000 

70.0 

.024 

.002 

0.000 

75.0 

0.000 

0.000 

0.000 

80.0 

0.000 

0.000 

0.000 

85.0 

0.000 

0.000 

0.000 

90.0 

0.000 

0.000 

0.000 

D6 


/ 


i 


i 


THERMAL  MODEL  INPUT  DATA  FOR  OAK  HICKORY 


GEOMETRICAL  VIEW 

ANGLE 

FACTOR  MATRIX 

ZENITH 

ANGLE 

5.0 

15.0 

25.0  35.0 

45.0 

55.0  65.0 

75.0 

85.0 

LAYER  1  1.0000 

.994? 

.9774  .9642 

.9573 

.9545  .9536 

.9536 

.  9536 

LAYER  2  0.0000 

.0019 

.0068  .0)05 

.0124 

.0131  .0134 

.0134 

.0134 

LAYER  3  0.0000 

.001 1 

.0047  .0075 

.0090 

.0096  .0098 

.0098 

.0098 

GROUND  0.0000 

.0023 

.0110  .0178 

.0213 

.0227  .0232 

.0232 

.0232 

LONG  WAVE 

TRANSFER  MATRIX 

FROM 

TO 

SKY 

LAYER  t  LAYER  2 

LAYER  3  GROUND 

LAYER  1 

.1595 

1.6741 

.0470 

.0338 

.0788 

LAYER  2 

.0281 

.7914 

.3539 

.2589 

.5607 

LAYER  3 

.0201 

.5442 

.2574 

.3496 

.3217 

AVERAGE  SHORTWAVE  ABSORPTION  COEFFICIENTS 


LAYER  1 

LAYER  2 

LAYER  3 

.089 

.042 

.040 

STOMATAL  RESISTANCE 
.0?  (MIN/CM) 


D7 


ENVIRONMENTAL  INPUT  DATA 


UALKER  BRANCH,  TENNESSEE  13  AUGUST  1979 
TIHE  AIR  TEMP  GRND  TEMP  UIND  SPEED  REL  HUM  GLOBAL  SUR 


(HOURS) 

(DEG  C) 

(DEG  C) 

100 

19.5 

19.5 

200 

19.3 

19.5 

300 

18.8 

19.4 

400 

18.3 

19.3 

500 

18.0 

19.2 

600 

17.7 

19.1 

700 

17.9 

19.1 

800 

19.3 

19.2 

900 

21.0 

19.3 

1000 

22.4 

19.5 

1100 

24.4 

19.9 

1200 

26.0 

20.4 

1300 

27.1 

20.7 

1400 

28.2 

21.0 

1500 

29.2 

21.4 

1600 

28.5 

21.6 

1700 

28.5 

21.6 

1800 

27.7 

21.7 

1900 

26.2 

21.6 

2000 

24.9 

21.4 

2100 

24.3 

21.3 

2200 

23.5 

21.2 

2300 

22.8 

21.0 

2400 

22.1 

20.7 

(M/SEC)  <U/M**2) 


3.1 

.83 

0.00 

3.2 

.85 

0.00 

2.3 

.87 

0.00 

2.3 

.91 

0.00 

2.7 

.94 

0.00 

2.4 

.96 

0.00 

3.0 

.96 

38.20 

3.0 

.93 

175.40 

2.4 

.88 

329.90 

4.3 

.84 

445.70 

3.5 

.77 

661.90 

3.6 

.71 

681.10 

2.8 

.65 

614.00 

2.9 

.64 

770.10 

2.8 

.61 

787.00 

2.8 

.61 

531.70 

2.7 

.62 

474.00 

2.4 

.65 

269.10 

2.3 

.70 

85.00 

2.2 

.75 

2.90 

2.9 

.76 

0.00 

3.0 

.78 

0.00 

3.2 

.80 

0.00 

3.4 

.83 

0.00 

D8 


/ 


ENVIRONMENTAL  INPUT  DATA 


UALKER  BRANCH,  TENNESSEE  19  AUGUST  1979 


TIME 

(HOURS) 

AIR  TEMP 
(DEG  C) 

GRND  TEMP 
(DEG  C) 

WIND  SPEED 
(M/SEC) 

REL  HUM 

100 

22.1 

20.7 

3.4 

.85 

200 

21.5 

20.6 

3.5 

.91 

300 

21.1 

20.6 

3.2 

.92 

400 

20.4 

20.5 

2.4 

.99 

500 

19.8 

20.3 

2.6 

1.00 

600 

19.6 

20.3 

1.8 

1.00 

700 

19.1 

20.2 

1.5 

1 .00 

800 

20.8 

20.3 

1.7 

1.00 

900 

23.8 

20.4 

1.5 

.91 

1000 

26.1 

20.6 

2.4 

.80 

1100 

27.8 

21.0 

2.5 

.73 

1200 

29.2 

21.5 

2.8 

.69 

1300 

30.4 

21.8 

2.2 

.63 

1400 

31.7 

22.2 

2.1 

.61 

1500 

31.3 

22.6 

2.1 

.60 

1600 

31.1 

22.8 

2.1 

.61 

1700 

30.2 

22.8 

2.5 

.63 

1800 

29.5 

22.8 

1.7 

.67 

1900 

27.9 

22.7 

1.9 

.74 

2000 

26.2 

22.4 

2.6 

.71 

2100 

25.9 

22.2 

2.1 

.69 

2200 

25.1 

22.0 

2.3 

.74 

2300 

24.3 

22.9 

2.0 

.81 

2400 

22.1 

21.7 

2.0 

.81 

GLOBAL  SUR 
<  U/M*<2) 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

22.80 

182.30 

365.00 

582.10 

751 .60 

753.60 
827.50 

917.40 
778.00 

620.40 

457.60 
251.70 

81.40 

2.60 

0.00 

0.00 

0.00 

0.00 


In  accordance  with  letter  from  DAEN-RDC,  DAEN-ASI  dated 
22  July  1977,  Subject:  Facsimile  Catalog  Cards  for 
Laboratory  Technical  Publications,  a  facsimile  catalog 
card  in  Library  of  Congress  MARC  format  is  reproduced 
below. 


Thermal  vegetation  canopy  model  studies  :  final  report  /  by 
J.A.  Smith  ...  [et  al ] ,  (Department  of  Wood  Science, 
College  of  Forestry  and  Natural  Resources,  Colorado 
State  University).  —  Vicksburg,  Miss.  :  U.S.  Army 
Engineer  Waterways  Experiment  Station  ;  Springfield, 

Va.  :  available  from  NTIS,  [I98l]. 

213  p.  in  various  pagings  :  ill.  ;  27  cm.  —  (Technical 
report  /  U.S.  Army  Engineer  Waterways  Experiment  Station  ; 
EL-81 -6) 

Cover  title. 

"August  1981." 

"Prepared  for  Headquarters,  Department  of  the  Army, 
under  Project  No.  l*A762730ATl*2,  Task  Al*,  Work  Unit  003 
(Contract  No.  DACW  39-77-C-0073) 

"Monitored  by  Environmental  Laboratory,  U.S.  Army 
Engineer  Waterways  Experiment  Station." 

Bibliography:  p.  !*1*-1*5. 


Thermal  vegetation  canopy  model  studies  :  ...  1981. 

(Card  2) 

1.  Computer  simulation.  2.  Infra-red  detectors. 

3.  Remote  sensing.  U.  Thermal  analysis.  5.  Vegetation 
classification.  I.  Smith,  J.A.  II.  Colorado  State 
University.  College  of  Forestry  and  Natural  Resources. 
III.  United  States.  Department  of  the  Army.  IV.  U.S. 

Army  Engineer  Waterways  Experiment  Station.  Environmental 
Laboratory.  IV.  Series:  Technical  report  (U.S.  Army 
Engineer  Waterways  Experiment  Station)  ;  EL-81-6. 

TA7.W3L  no. EL-81-6 


